Pattern of litterfall production throughout a tropical humid forest chronosequence
Judith Karine Obando Guzmán, Luis C. Beltrán*, Carlos Hernando Rodríguez Lilia L. Roa-Fuentes
2023-06-19
Setup
Loading Packages
We will start by loading some useful packages. Using the
library function, we can load these packages if they have
previously been installed in your R. If they have not been
previously installed, replace the library function with
install.packages function and place the package name in
quotation marks. That will install the package. Next, run the
library function again and it will load the package. You
only need to install a package once.
library(plyr) # Necessary for Rmisc
library(Rmisc) # Necessary for summarySE
library(dplyr) # Necessary for group_by
library(reshape) # Necessary for melt & cast
library(ggplot2) # Necessary for ggplots
library(grid) # Necessary for pushViewPort
library(ggpubr) # Necessary for ggarrange
library(conover.test) # Necessary for conover test
library(lmtest) # Necessary for homoscedasticity test.
library(lme4) # Generalized linear models and mixed effects.
library(lmerTest) # Necessary to get p-vals in glms
library(lubridate) # Necessary for monthiness?
library(ggthemes) # Necessary for cool themes
library(emmeans) # Display all pairwise comparisons
library(ggfortify) # Necessary for PCA.
library(patchwork) # Necessary for mixing ggplots and legends
library(kableExtra) # Necessary for 'kbl' function.
library(report) # Necessary for 'sessionInfo' function.Reading Dataframes
We will set the working directory using the setwd
function. If you are looking to run this code in your computer, be sure
to change the path to your data. load dataframes using the
read.csv function.
setwd("~/Desktop/R")
LeafLitt <- read.csv("Data/Leaf.Litter.Ecosystems.csv")
Clim <- read.csv("Data/Clim.csv")
Struc <- read.csv("Data/StrucDF.csv")
Solar <- read.csv("Data/SolarRadiation.csv")
Irad <- read.csv("Data/SolarIrradiation.csv")
Composition <- read.csv("Data/PlotStructureCaqueta.csv") Now, the next step is to get these dataframes in the format that we
want. Several columns are of the character class but we want them as
factors. Here the across function selects the columns to be
converted to factors based on their names, and as.factor is
applied to each selected column using the mutate function
from the dplyr package. Finally, the resulting factors are
assigned back to the corresponding columns using the %>%
pipe operator.
Correcting Variable Class
LeafLitt <- LeafLitt %>%
dplyr::mutate(across(c("Paisaje", "Sample", "Locality", "Successional.State", "Succession",
"Plot", "Record"), as.factor))
Clim <- Clim %>%
dplyr::mutate(across(c("Landscape", "Plot", "Sample"), as.factor))
Struc <- Struc %>%
dplyr::mutate(across(c("Landscape", "Plot", "Succession", "CodeZ"), as.factor))
Solar <- Solar %>%
dplyr::mutate(across(c("Estacion", "Municipio"), as.factor))
Irad$Month <- as.factor(Irad$Month) # Florencia Data.
Composition <- Composition %>%
dplyr::mutate(across(c("Landscape", "Chrono", "Chronosecuence", "Plot", "Species"), as.factor))Great. Next we need to change some values from Muestro 5 (353 (DM2), 365(RM), and 377(RM)) into NAs, given that Muestreo 5 is incomplete for some plots.
Setting NAs
LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA
LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "377")] <- NADataframe Setup
In the LeafLitt dataframe, there is one level of the
variable successiont hat we will be removing, DV1. We will
create a new dataframe that excludes this level using the
droplevels and subset functions.
LeafLitt2 <- droplevels(subset(LeafLitt, LeafLitt$Succession != "DV1"))Next, we will use the paste function to create a new
variable that is a combination of the text from the
Paisaje, Successional.State, and
Plot (CodeX). We will also create one for a
combination of Succession and Plot
(CodeZ). These variables will be used for matching
dataframes further down the code.
LeafLitt2$CodeX <- paste(LeafLitt2$Paisaje , LeafLitt2$Successional.State, LeafLitt2$Plot)
LeafLitt2$CodeZ <- paste(LeafLitt2$Succession , LeafLitt2$Plot)Here we will create a CodeX for the Clim
dataframe as well. Then we will use CodeX to attach the
Succession column from LeafLitt2 using the
with and match functions. Next we remove NA
values from Clim dataframe using droplevels
and subset functions to create a new dataframe,
Clim2. A CodeY variable is created by pasting
CodeX and Sample for Clim2 and
LeafLitt2.
Clim$CodeX <- paste(Clim$Landscape , Clim$Successional.State, Clim$Plot)
Clim$Succession <- with(LeafLitt2, Succession[match(Clim$CodeX, CodeX)])
Clim2 <- droplevels(subset(Clim, Clim$Succession != "NA"))
Clim2$CodeY <- paste(Clim2$CodeX , Clim2$Sample)
LeafLitt2$CodeY <- paste(LeafLitt2$CodeX , LeafLitt2$Sample)Using the shared CodeY, we will now insert the
precipitation and temperature data into the LeafLitter2
dataframe. Note that Rain1 has a 1-mth lag,
Rain2 has 2-mth lag, and Rain3 3-mth long.
Same for temperature.
LeafLitt2$RainNow <- with(Clim2, RainNow[match(LeafLitt2$CodeY, Clim2$CodeY)])
LeafLitt2$Rain1 <- with(Clim2, Rain1[match(LeafLitt2$CodeY, Clim2$CodeY)])
LeafLitt2$Rain2 <- with(Clim2, Rain2[match(LeafLitt2$CodeY, Clim2$CodeY)])
LeafLitt2$Rain3 <- with(Clim2, Rain3[match(LeafLitt2$CodeY, Clim2$CodeY)])
LeafLitt2$TempNow <- with(Clim2, TempNow[match(LeafLitt2$CodeY, Clim2$CodeY)])
LeafLitt2$Temp1 <- with(Clim2, Temp1[match(LeafLitt2$CodeY, Clim2$CodeY)])
LeafLitt2$Temp2 <- with(Clim2, Temp2[match(LeafLitt2$CodeY, Clim2$CodeY)])
LeafLitt2$Temp3 <- with(Clim2, Temp3[match(LeafLitt2$CodeY, Clim2$CodeY)]) We’ll use the same trick as before to make our Code variables into factors.
LeafLitt2 <- LeafLitt2 %>%
dplyr::mutate(across(c("CodeX", "CodeZ", "CodeY"), as.factor))For various analyses we will want to deal with landscapes separately,
so here we create one dataframe for Hill plots (Lomerio)
and for Mountain using droplevels and
subset functions.
Lomerio <- droplevels(subset(LeafLitt2, LeafLitt2$Paisaje == "Lomerio"))
Mountain <- droplevels(subset(LeafLitt2, LeafLitt2$Paisaje == "Mountain"))Now, the months in these dataframes are coded as numbers, and these numbers differ by dataframe, so for each dataframe, we will manually relabel the months.
levels(Lomerio$Sample) <- list("Aug" = "1", "Sept" = "2", "Oct" = "3",
"Nov" = "4", "Dec" = "5", "Jan" = "6",
"Feb" = "7", "Mar" = "8", "Apr" = "9",
"May" = "10", "Jun" = "11", "Jul" = "12")
levels(Mountain$Sample) <- list("Sept" = "1", "Oct" = "2", "Nov" = "3",
"Dec" = "4", "Jan" = "5", "Feb" = "6",
"Mar" = "7", "Apr" = "8", "May" = "9",
"Jun" = "10", "Jul" = "11", "Aug" = "12")We continue splitting, this time by Age Group.
LomeDL1 <- droplevels(subset(Lomerio, Lomerio$Succession == "DL1"))
LomeDL2 <- droplevels(subset(Lomerio, Lomerio$Succession == "DL2"))
LomeRL <- droplevels(subset(Lomerio, Lomerio$Succession == "RL"))
MounDM1 <- droplevels(subset(Mountain, Mountain$Succession == "DM1"))
MounDM2 <- droplevels(subset(Mountain, Mountain$Succession == "DM2"))
MounRL <- droplevels(subset(Mountain, Mountain$Succession == "RM"))Descriptive Statistics
Mean and Standard Error
We can obtain mean and standard error for each variable in each
landscape by running the summarySE function on the
LeafLitt2 dataframe and setting Paisaje as the
grouping variable. Before that, we will use the options
function to set how many digits we want to see. We will wrap this up
within the kbl function, which allows me to present the
table in a neater format. We will also exclude some other columns using
the square brackets.
options(digits=3)
kbl(summarySE(LeafLitt2, measurevar="LeafLitter.Mg.ha.month", groupvars=c("Paisaje"), na.rm = TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Paisaje | N | LeafLitter.Mg.ha.month | se |
|---|---|---|---|
| Lomerio | 156 | 0.603 | 0.022 |
| Mountain | 213 | 0.433 | 0.015 |
kbl(summarySE(LeafLitt2, measurevar="Branches.Mg.ha.month", groupvars=c("Paisaje"), na.rm = TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Paisaje | N | Branches.Mg.ha.month | se |
|---|---|---|---|
| Lomerio | 156 | 0.040 | 0.003 |
| Mountain | 213 | 0.051 | 0.003 |
kbl(summarySE(LeafLitt2, measurevar="Detritus.Mg.ha.month", groupvars=c("Paisaje"), na.rm = TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Paisaje | N | Detritus.Mg.ha.month | se |
|---|---|---|---|
| Lomerio | 156 | 0.132 | 0.007 |
| Mountain | 213 | 0.135 | 0.006 |
kbl(summarySE(LeafLitt2, measurevar="Flowers.Mg.ha.month", groupvars=c("Paisaje"), na.rm = TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Paisaje | N | Flowers.Mg.ha.month | se |
|---|---|---|---|
| Lomerio | 156 | 0.002 | 0.000 |
| Mountain | 213 | 0.003 | 0.001 |
kbl(summarySE(LeafLitt2, measurevar="Fruits.Mg.ha.month", groupvars=c("Paisaje"), na.rm =
TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Paisaje | N | Fruits.Mg.ha.month | se |
|---|---|---|---|
| Lomerio | 156 | 0.022 | 0.003 |
| Mountain | 213 | 0.018 | 0.002 |
kbl(summarySE(LeafLitt2, measurevar="Total.Mg.ha.month", groupvars=c("Paisaje"), na.rm =
TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Paisaje | N | Total.Mg.ha.month | se |
|---|---|---|---|
| Lomerio | 156 | 0.807 | 0.028 |
| Mountain | 213 | 0.647 | 0.021 |
Percent Composition
Ok, so how much of the capture leaf litter is made of Leaf, Branch, Fruit, Flower, or Detritus?
options(digits=4)
DFX1 <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month) ~ Paisaje, data = LeafLitt2, FUN = sum, na.action=na.exclude)
DFX1$LeafP <- DFX1[,2]/DFX1[,7]*100
DFX1$BranchP <- DFX1[,3]/DFX1[,7]*100
DFX1$FlowerP <- DFX1[,4]/DFX1[,7]*100
DFX1$FruitP <- DFX1[,5]/DFX1[,7]*100
DFX1$DetritusP <- DFX1[,6]/DFX1[,7]*100
kbl(DFX1[c(1,8:12)], format = "html", table.attr = "style='width:50%;'") %>% kable_classic()| Paisaje | LeafP | BranchP | FlowerP | FruitP | DetritusP |
|---|---|---|---|---|---|
| Lomerio | 74.73 | 5.007 | 0.2344 | 2.717 | 16.38 |
| Mountain | 66.98 | 7.920 | 0.5166 | 2.800 | 20.93 |
DFX4 <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month) ~ Succession,
data = Lomerio, FUN = sum, na.action=na.exclude)
DFX4$LeafP <- DFX4[,2]/DFX4[,7]*100
DFX4$BranchP <- DFX4[,3]/DFX4[,7]*100
DFX4$FlowerP <- DFX4[,4]/DFX4[,7]*100
DFX4$FruitP <- DFX4[,5]/DFX4[,7]*100
DFX4$DetritusP <- DFX4[,6]/DFX4[,7]*100
kbl(DFX4[c(1,8:12)], format = "html", table.attr = "style='width:50%;'") %>% kable_classic()| Succession | LeafP | BranchP | FlowerP | FruitP | DetritusP |
|---|---|---|---|---|---|
| DL1 | 77.01 | 3.654 | 0.0208 | 3.342 | 14.99 |
| DL2 | 74.04 | 5.152 | 0.2038 | 3.072 | 16.61 |
| RL | 74.66 | 5.376 | 0.3598 | 2.017 | 16.66 |
DFX4 <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month) ~ Succession,
data = Mountain, FUN = sum, na.action=na.exclude)
DFX4$LeafP <- DFX4[,2]/DFX4[,7]*100
DFX4$BranchP <- DFX4[,3]/DFX4[,7]*100
DFX4$FlowerP <- DFX4[,4]/DFX4[,7]*100
DFX4$FruitP <- DFX4[,5]/DFX4[,7]*100
DFX4$DetritusP <- DFX4[,6]/DFX4[,7]*100
kbl(DFX4[c(1,8:12)], format = "html", table.attr = "style='width:50%;'") %>% kable_classic()| Succession | LeafP | BranchP | FlowerP | FruitP | DetritusP |
|---|---|---|---|---|---|
| DM1 | 71.71 | 6.694 | 0.3250 | 1.651 | 18.70 |
| DM2 | 64.04 | 8.865 | 0.4382 | 2.915 | 22.93 |
| RM | 66.74 | 7.797 | 0.7298 | 3.477 | 20.41 |
Diagnostic Statistics
Normality Test
We will use Shapiro-Wilk test, available through the
shapiro.test function, to test for normality in our
response variables. The Shapiro-Wilk test generates a test statistic,
usually denoted as W. The closer the value of W is to 1, the more the
data resemble a normal distribution. The shapiro.test
function also provides a p-value, which represents the probability of
obtaining the observed data if the null hypothesis (data are normally
distributed) is true. If the p-value is greater than the chosen
significance level (commonly 0.05), we fail to reject the null
hypothesis, suggesting that the data can be considered approximately
normally distributed. If the p-value is less than the chosen
significance level, we reject the null hypothesis and conclude that the
data significantly deviate from a normal distribution.
shapiro.test(LeafLitt2$LeafLitter.Mg.ha.month) ##
## Shapiro-Wilk normality test
##
## data: LeafLitt2$LeafLitter.Mg.ha.month
## W = 0.94, p-value = 5e-11shapiro.test(LeafLitt2$Branches.Mg.ha.month) ##
## Shapiro-Wilk normality test
##
## data: LeafLitt2$Branches.Mg.ha.month
## W = 0.85, p-value <2e-16shapiro.test(LeafLitt2$Detritus.Mg.ha.month) ##
## Shapiro-Wilk normality test
##
## data: LeafLitt2$Detritus.Mg.ha.month
## W = 0.85, p-value <2e-16shapiro.test(LeafLitt2$Flowers.Mg.ha.month) ##
## Shapiro-Wilk normality test
##
## data: LeafLitt2$Flowers.Mg.ha.month
## W = 0.42, p-value <2e-16shapiro.test(LeafLitt2$Fruits.Mg.ha.month) ##
## Shapiro-Wilk normality test
##
## data: LeafLitt2$Fruits.Mg.ha.month
## W = 0.61, p-value <2e-16shapiro.test(LeafLitt2$Total.Mg.ha.month) ##
## Shapiro-Wilk normality test
##
## data: LeafLitt2$Total.Mg.ha.month
## W = 0.95, p-value = 2e-09None of the variables follow a normal distribution. This isn’t necessarily wrong or bad, but it does mean that we may obtain more accurate models if we account for lack of normality.
Data Distributions
Let’s take a look at some distributions. For this we will create a
histogram and kernel density estimate, which is basically a smoothed out
histogram using the geom_histogram and
geom_density arguments of ggplot. We will also
use geom_vline to create a vertical line where the mean of
our species richness is.
ggplot(LeafLitt2, aes(x=LeafLitter.Mg.ha.month)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(LeafLitter.Mg.ha.month, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Leaf Litter (Mg/Ha/Mth") +
ylab("Density") +
theme_bw(base_size = 15)
ggplot(LeafLitt2, aes(x=Branches.Mg.ha.month)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Branches.Mg.ha.month, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Branch Litter (Mg/Ha/Mth") +
ylab("Density") +
theme_bw(base_size = 15)
ggplot(LeafLitt2, aes(x=Detritus.Mg.ha.month)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Detritus.Mg.ha.month, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Detritus Litter (Mg/Ha/Mth") +
ylab("Density") +
theme_bw(base_size = 15)
ggplot(LeafLitt2, aes(x=Flowers.Mg.ha.month)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Flowers.Mg.ha.month, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Flower Litter (Mg/Ha/Mth") +
ylab("Density") +
theme_bw(base_size = 15)
ggplot(LeafLitt2, aes(x=Fruits.Mg.ha.month)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Fruits.Mg.ha.month, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Fruit Litter (Mg/Ha/Mth") +
ylab("Density") +
theme_bw(base_size = 15)
ggplot(LeafLitt2, aes(x=Total.Mg.ha.month)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Total.Mg.ha.month, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Total Litter (Mg/Ha/Mth") +
ylab("Density") +
theme_bw(base_size = 15)
All distributions are right-skewed, some more than others.
Homoscedasticity & Linearity
Next let’s check for homoscedasticity and linearity. To do this, we’l create our first linear models.
lm1 <- lm(LeafLitter.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm1, 1) # Residual vs Fitted looks fine.
plot(lm1, 3) # Scale-Location looks fine too.
lm2 <- lm(LeafLitter.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm2, 1) # Residual vs Fitted looks fine.
plot(lm2, 3) # Scale-Location looks fine too.
lm3 <- lm(Branches.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm3, 1) # Residual vs Fitted looks fine.
plot(lm3, 3) # Scale-Location shows increase of residuals with fitted values.
lm4 <- lm(Branches.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm4, 1) # Residual vs Fitted looks fine.
plot(lm4, 3) # Scale-Location looks fine too.
lm5 <- lm(Flowers.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm5, 1) # Residual vs Fitted looks fine.
plot(lm5, 3) # Scale-Location Shows mild increase.
lm6 <- lm(Flowers.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm6, 1) # Residual vs Fitted looks fine.
plot(lm6, 3) # Scale-Location shows increase of residuals with fitted values.
lm7 <- lm(Fruits.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm7, 1) # Residual vs Fitted looks fine.
plot(lm7, 3) # Scale-Location looks fine too.
lm8 <- lm(Fruits.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm8, 1) # Residual vs Fitted looks fine.
plot(lm8, 3) # Scale-Location looks fine too.
lm9 <- lm(Detritus.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm9, 1) # Residual vs Fitted looks fine.
plot(lm9, 3) # Scale-Location looks fine too.
lm10 <- lm(Detritus.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm10, 1) # Residual vs Fitted looks fine.
plot(lm10, 3) # Scale-Location looks fine too.
lm11 <- lm(Total.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm11, 1) # Residual vs Fitted looks fine.
plot(lm11, 3) # Scale-Location looks fine too.
lm12 <- lm(Total.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm12, 1) # Residual vs Fitted looks fine.
plot(lm12, 3) # Scale-Location looks fine too. 
Scedasticity and linearity are largely fine. We seem to mostly have a lack of non-normality with a right-skew.
Transformations
For left-skewed data—tail is on the left, negative skew, common transformations include (i) square root, (ii) cube root, and (iii) log.
First test: Square root transformation.
hist(sqrt(LeafLitt2$LeafLitter.Mg.ha.month))
hist(sqrt(LeafLitt2$Branches.Mg.ha.month))
hist(sqrt(LeafLitt2$Detritus.Mg.ha.month))
hist(sqrt(LeafLitt2$Flowers.Mg.ha.month))
hist(sqrt(LeafLitt2$Fruits.Mg.ha.month))
hist(sqrt(LeafLitt2$Total.Mg.ha.month))
Square-root seems to work well for Leaf, Detritus, and Total. Branches
sorta.
Second test: Log transformation.
hist(log(LeafLitt2$Flowers.Mg.ha.month)) # Great
hist(log(LeafLitt2$Fruits.Mg.ha.month)) # Ok
hist(log(LeafLitt2$Branches.Mg.ha.month)) # Right skew.
Log works well for flowers and fruit. Branches needs something else.
hist(LeafLitt2$Branches.Mg.ha.month^(1/3)) 
Cube root best for Branches.
Note: Using the models without transformations does not yield radically different conclusions. However, we maintain the use of transformations to remain consistent with the literature and with reccommendations from biostatisticians.
Descriptive Results
Flowering Peaks
Flower litter, when does it fall with greater abundance? Let’s find
out for each age group in each landscape. First for Hills. Let’s
calculate means of Flower production by age Month using the
aggregate function wrapped within a kbl
function for style.
options(digits=7)
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = LomeRL, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()| Sample | Flowers.Mg.ha.month |
|---|---|
| Aug | 0.02000 |
| Sept | 0.00150 |
| Oct | 0.00250 |
| Nov | 0.01000 |
| Dec | 0.00025 |
| Jan | 0.00025 |
| Feb | 0.00100 |
| Mar | 0.00050 |
| Apr | 0.00075 |
| May | 0.00025 |
| Jun | 0.00200 |
| Jul | 0.00350 |
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = LomeDL2, FUN = mean, na.rm = TRUE), format = "html",
table.attr = "style='width:30%;'") %>% kable_classic() | Sample | Flowers.Mg.ha.month |
|---|---|
| Aug | 0.0018333 |
| Sept | 0.0006667 |
| Oct | 0.0015000 |
| Nov | 0.0013333 |
| Dec | 0.0025000 |
| Jan | 0.0010000 |
| Feb | 0.0055000 |
| Mar | 0.0006667 |
| Apr | 0.0028333 |
| May | 0.0001667 |
| Jun | 0.0000000 |
| Jul | 0.0021667 |
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = LomeDL1, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:30%;'") %>% kable_classic() | Sample | Flowers.Mg.ha.month |
|---|---|
| Aug | 0.0000000 |
| Sept | 0.0000000 |
| Oct | 0.0000000 |
| Nov | 0.0000000 |
| Dec | 0.0003333 |
| Jan | 0.0000000 |
| Feb | 0.0006667 |
| Mar | 0.0000000 |
| Apr | 0.0000000 |
| May | 0.0003333 |
| Jun | 0.0000000 |
| Jul | 0.0000000 |
ggplot(data=Lomerio, aes(x=Sample, y=Flowers.Mg.ha.month )) + # Flowers
geom_bar(stat="identity", aes(fill = Succession), position="dodge") +
xlab("Month") +
ylab(bquote('Flower Litter '(Mg~ha^-1~m^-1))) +
scale_fill_manual(values=c("#82C27B", "#43833C", "#23431F")) +
guides(fill=guide_legend(title="Age Group")) +
theme_bw(base_size = 12)
Next for the mountain landscape.
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = MounRL, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic() | Sample | Flowers.Mg.ha.month |
|---|---|
| Sept | 0.0146667 |
| Oct | 0.0098333 |
| Nov | 0.0020000 |
| Dec | 0.0010000 |
| Jan | 0.0037500 |
| Feb | 0.0000000 |
| Mar | 0.0028333 |
| Apr | 0.0010000 |
| May | 0.0008333 |
| Jun | 0.0036667 |
| Jul | 0.0051667 |
| Aug | 0.0181667 |
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = MounDM2, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic() | Sample | Flowers.Mg.ha.month |
|---|---|
| Sept | 0.0026667 |
| Oct | 0.0035000 |
| Nov | 0.0036667 |
| Dec | 0.0035000 |
| Jan | 0.0010000 |
| Feb | 0.0018333 |
| Mar | 0.0038333 |
| Apr | 0.0005000 |
| May | 0.0075000 |
| Jun | 0.0063333 |
| Jul | 0.0023333 |
| Aug | 0.0015000 |
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = MounDM1, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic() | Sample | Flowers.Mg.ha.month |
|---|---|
| Sept | 0.0010000 |
| Oct | 0.0001667 |
| Nov | 0.0003333 |
| Dec | 0.0018333 |
| Jan | 0.0091667 |
| Feb | 0.0008333 |
| Mar | 0.0016667 |
| Apr | 0.0000000 |
| May | 0.0005000 |
| Jun | 0.0005000 |
| Jul | 0.0021667 |
| Aug | 0.0008333 |
ggplot(data=Mountain, aes(x=Sample, y=Flowers.Mg.ha.month )) + # Flowers
geom_bar(stat="identity", aes(fill = Succession), position="dodge") +
xlab("Month") +
ylab(bquote('Flower Litter '(Mg~ha^-1~m^-1))) +
scale_fill_manual(values=c("#FBBA84", "#CE803F", "#67360D")) +
guides(fill=guide_legend(title="Age Group")) +
theme_bw(base_size = 10)
Fruiting Peaks
options(digits=5)
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = LomeRL, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()| Sample | Fruits.Mg.ha.month |
|---|---|
| Aug | 0.00400 |
| Sept | 0.02325 |
| Oct | 0.00925 |
| Nov | 0.02200 |
| Dec | 0.01150 |
| Jan | 0.02550 |
| Feb | 0.02700 |
| Mar | 0.03500 |
| Apr | 0.01650 |
| May | 0.01925 |
| Jun | 0.02450 |
| Jul | 0.02050 |
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = LomeDL2, FUN = mean, na.rm = TRUE), format = "html",
table.attr = "style='width:30%;'") %>% kable_classic() | Sample | Fruits.Mg.ha.month |
|---|---|
| Aug | 0.00783 |
| Sept | 0.05217 |
| Oct | 0.02217 |
| Nov | 0.00967 |
| Dec | 0.01000 |
| Jan | 0.05083 |
| Feb | 0.03567 |
| Mar | 0.01717 |
| Apr | 0.02867 |
| May | 0.01883 |
| Jun | 0.00933 |
| Jul | 0.04167 |
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = LomeDL1, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:30%;'") %>% kable_classic() | Sample | Fruits.Mg.ha.month |
|---|---|
| Aug | 0.00300 |
| Sept | 0.00167 |
| Oct | 0.00500 |
| Nov | 0.00333 |
| Dec | 0.01367 |
| Jan | 0.00667 |
| Feb | 0.00467 |
| Mar | 0.01433 |
| Apr | 0.02933 |
| May | 0.04300 |
| Jun | 0.00767 |
| Jul | 0.08167 |
ggplot(data=Lomerio, aes(x=Sample, y=Fruits.Mg.ha.month )) + # Fruits
geom_bar(stat="identity", aes(fill = Succession), position="dodge") +
xlab("Month") +
ylab(bquote('Fruit Litter '(Mg~ha^-1~m^-1))) +
scale_fill_manual(values=c("#82C27B", "#43833C", "#23431F")) +
guides(fill=guide_legend(title="Age Group")) +
theme_bw(base_size = 12)
Next for the mountain landscape.
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = MounRL, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic() | Sample | Fruits.Mg.ha.month |
|---|---|
| Sept | 0.04417 |
| Oct | 0.04500 |
| Nov | 0.02167 |
| Dec | 0.01583 |
| Jan | 0.01225 |
| Feb | 0.01567 |
| Mar | 0.01833 |
| Apr | 0.02017 |
| May | 0.02583 |
| Jun | 0.01867 |
| Jul | 0.01383 |
| Aug | 0.04650 |
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = MounDM2, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic() | Sample | Fruits.Mg.ha.month |
|---|---|
| Sept | 0.06767 |
| Oct | 0.02100 |
| Nov | 0.01017 |
| Dec | 0.01200 |
| Jan | 0.01080 |
| Feb | 0.01483 |
| Mar | 0.01000 |
| Apr | 0.01317 |
| May | 0.02917 |
| Jun | 0.02350 |
| Jul | 0.01333 |
| Aug | 0.02900 |
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = MounDM1, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic() | Sample | Fruits.Mg.ha.month |
|---|---|
| Sept | 0.00683 |
| Oct | 0.00283 |
| Nov | 0.00333 |
| Dec | 0.00133 |
| Jan | 0.00617 |
| Feb | 0.01167 |
| Mar | 0.01050 |
| Apr | 0.00567 |
| May | 0.00733 |
| Jun | 0.00933 |
| Jul | 0.00700 |
| Aug | 0.02450 |
ggplot(data=Mountain, aes(x=Sample, y=Fruits.Mg.ha.month )) + # Fruits
geom_bar(stat="identity", aes(fill = Succession), position="dodge") +
xlab("Month") +
ylab(bquote('Fruit Litter '(Mg~ha^-1~m^-1))) +
scale_fill_manual(values=c("#FBBA84", "#CE803F", "#67360D")) +
guides(fill=guide_legend(title="Age Group")) +
theme_bw(base_size = 10)
Reference Age Groups
Lomerio reference plots appear to have higher leaf litter and
structural variable totals than Mountain reference plots. Let’s check.
First we’ll create a new dataframe that holds only reference plots and
we’ll call it ReferencePlots. Then, we will use the
aggregate function to obtain means for each variable of
interest, setting Age Groups (Succession) as the grouping
variable.
options(digits=5)
ReferencePlots <- droplevels(subset(LeafLitt2, LeafLitt2$Succession == c("RM", "RL")))
kbl(aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month) ~ Succession,data = ReferencePlots, FUN = mean, na.action=na.exclude), format = "html", table.attr = "style='width:100%;'", row.names = FALSE) %>%
kable_classic()| Succession | LeafLitter.Mg.ha.month | Branches.Mg.ha.month | Flowers.Mg.ha.month | Fruits.Mg.ha.month | Detritus.Mg.ha.month | Total.Mg.ha.month |
|---|---|---|---|---|---|---|
| RL | 0.75613 | 0.05987 | 0.00267 | 0.02425 | 0.18196 | 1.03412 |
| RM | 0.46803 | 0.05847 | 0.00494 | 0.02329 | 0.14229 | 0.70279 |
ReferenceStrucPlots <- droplevels(subset(Struc, Struc$Succession == c("RM", "RL")))
kbl(aggregate(cbind(TreesHA, BasalHA, Biomass, Richness, Simpson_1.D) ~ Landscape,
data = ReferenceStrucPlots, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:100%;'", row.names = FALSE) %>% kable_classic()| Landscape | TreesHA | BasalHA | Biomass | Richness | Simpson_1.D |
|---|---|---|---|---|---|
| Lomerio | 614.00 | 25.35 | 252.05 | 87.500 | 0.9755 |
| Mountain | 438.67 | 13.00 | 117.67 | 43.667 | 0.9300 |
Sure enough, Lomerio (RL) has higher leaf, detritus, and total leaf litter than mountain (RM). RL also has higher tree density, basal area, biomass, richness, and a lil bit of diversity.
Rainfall Differences
How do plots differ in rainfall? We can take a look using the same combination of functions as before. First we’ll take a look at the Hill landscape plots.
kbl(aggregate(cbind(RainNow) ~ Plot, data = Lomerio, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:20%;'", row.names = FALSE) %>% kable_classic()| Plot | RainNow |
|---|---|
| P10 | 308.03 |
| P11 | 302.92 |
| P12 | 308.61 |
| P13 | 308.61 |
| P14 | 324.03 |
| P16 | 323.96 |
| P2 | 313.99 |
| P3 | 308.44 |
| P4 | 296.90 |
| P5 | 302.92 |
| P6 | 313.42 |
Ok, in the Hillplots, rainfall dances around 300, a couple have 320 mm. Now for the mountain.
kbl(aggregate(cbind(RainNow) ~ Plot, data = Mountain, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:20%;'", row.names = FALSE) %>% kable_classic()| Plot | RainNow |
|---|---|
| P1 | 256.25 |
| P10 | 286.60 |
| P11 | 288.45 |
| P12 | 284.69 |
| P13 | 280.93 |
| P14 | 284.75 |
| P15 | 265.91 |
| P16 | 256.25 |
| P17 | 265.91 |
| P18 | 255.08 |
| P2 | 256.25 |
| P6 | 256.25 |
| P7 | 256.25 |
| P8 | 265.91 |
| P9 | 265.91 |
In the mountain landscape, precipitation (mm) values dance around 260, a couple have 288.
So now, do the landscapes differ in precipitation? Let’s find out via
anova.
anova(lm(data=LeafLitt2, RainNow ~ Paisaje))## Analysis of Variance Table
##
## Response: RainNow
## Df Sum Sq Mean Sq F value Pr(>F)
## Paisaje 1 140235 140235 9.33 0.0024 **
## Residuals 370 5561807 15032
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Indeed, landscapes signficantly differ in the amount of rainfall they receive. Hill landscape likely receives more, let’s see.
kbl(summarySE(LeafLitt2, measurevar="RainNow", groupvars=c("Paisaje"), na.rm = TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:40%;'", row.names = FALSE) %>% kable_classic()| Paisaje | N | RainNow | se |
|---|---|---|---|
| Lomerio | 156 | 310.33 | 10.7280 |
| Mountain | 216 | 270.98 | 7.7356 |
Sunlight and Irradiation
Before anything, we need to get our Sunlight and Irradiation
dataframes into formats that are easier to use. For this, we apply the
melt function to create SolarMelt, a dataframe
where we have a single column for Month, a single column
for Municipio, and a single column for
Sunlight.
SolarMelt <- melt(Solar,id.vars=c("Estacion", "Municipio"))
colnames(SolarMelt) <- c("Estacion", "Municipio", "Month", "Sunlight")Now let’s visualize using ggplot. Do different Municipios in Caqueta receive wildly different hours of sunlight? Likely no.
ggplot(SolarMelt, aes(x=Month, y=Sunlight, group = Estacion, color = Estacion)) +
geom_line(size=1, alpha=0.9, linetype=1) +
ylab("Sunlight (hrs/day)") +
scale_colour_manual(values = c("#00AFBB", "#E7B800", "#FC4E07","#FFDB6D", "#C4961A", "#F4EDCA",
"#D16103", "#C3D7A4", "#52854C", "#4E84C4", "#293352")) +
theme_bw(base_size = 12) 
Great. It appears that all municipios share the same general sunlight
pattern, so maybe we can simply use an average! Let’s extract the
average sunlight per station and save it as
SolarMeltAve.
SolarMeltAve <- droplevels(subset(SolarMelt, SolarMelt$Estacion == "Average"))Next, we need to do some setup to create a dataframe that holds
litter production and our climatic variables. First, we will use the
select function from the dplyr package to
specify the columns we wish to keep in our new dataframes,
Lome3 and Moun3. Next, we want to calculate
the average litter production and precipitation by month. These will
then be Lome3Ave and Moun3Ave. After renaming
columns with colnames, we merge with our sunlight dataframe
using the merge function, creating the AveLom
and AveMon dataframes.
Lome3 <- Lomerio %>% dplyr::select(Sample, LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Detritus.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Total.Mg.ha.month, Precipitation)
Moun3 <- Mountain %>% dplyr::select(Sample, LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Detritus.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Total.Mg.ha.month, Precipitation)
Lome3Ave <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month, Precipitation) ~ Sample,
data = Lome3, FUN = mean, na.action=na.exclude)
Moun3Ave <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month, Precipitation) ~ Sample,
data = Moun3, FUN = mean, na.action=na.exclude)
colnames(Lome3Ave) <- c("Month", "LeafLitter.Mg.ha.month", "Branches.Mg.ha.month", "Detritus.Mg.ha.month",
"Flowers.Mg.ha.month", "Fruits.Mg.ha.month", "Total.Mg.ha.month", "Precipitation")
colnames(Moun3Ave) <- c("Month", "LeafLitter.Mg.ha.month", "Branches.Mg.ha.month", "Detritus.Mg.ha.month",
"Flowers.Mg.ha.month", "Fruits.Mg.ha.month", "Total.Mg.ha.month", "Precipitation")
AveLom <- merge(SolarMeltAve, Lome3Ave, by="Month")
AveMon <- merge(SolarMeltAve, Moun3Ave, by="Month")Next, we merge these dataframes with the Irad dataframe
also through the merge function.
AveLomX <- merge(AveLom, Irad, by="Month")
AveMonX <- merge(AveMon, Irad, by="Month")Now lets quickly visualize how precipitation and sunlight behave within the hill and mountain landscapes.
ggplot(AveLom, aes(x=Precipitation, y=Sunlight)) +
geom_point(size=1, alpha=0.9, linetype=1) +
geom_smooth(method=lm, colour = "black", fill="#23431F") +
ylab("Sunlight (hrs/day)") +
theme_bw(base_size = 15)
ggplot(AveMon, aes(x=Precipitation, y=Sunlight)) +
geom_point(size=1, alpha=0.9, linetype=1) +
geom_smooth(method=lm, colour = "black", fill="#67360D") +
ylab("Sunlight (hrs/day)") +
theme_bw(base_size = 15) 
Good! Certainly look negatively correlated, but lets check using the
cor.test function. We’ll run both, Pearson’s and Spearman’s
Rank Correlation to check. We’ll also run the correaltions between
precipitation and Irradiation.
cor.test(AveLomX$Precipitation, AveLomX$Sunlight, method = "pearson")##
## Pearson's product-moment correlation
##
## data: AveLomX$Precipitation and AveLomX$Sunlight
## t = -2.57, df = 10, p-value = 0.028
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.884504 -0.089213
## sample estimates:
## cor
## -0.63082cor.test(AveLomX$Precipitation, AveLomX$Sunlight, method = "spearman")##
## Spearman's rank correlation rho
##
## data: AveLomX$Precipitation and AveLomX$Sunlight
## S = 460, p-value = 0.036
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.60702cor.test(AveMonX$Precipitation, AveMonX$Sunlight, method = "pearson")##
## Pearson's product-moment correlation
##
## data: AveMonX$Precipitation and AveMonX$Sunlight
## t = -2.2, df = 10, p-value = 0.052
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.8624520 0.0038045
## sample estimates:
## cor
## -0.57134cor.test(AveMonX$Precipitation, AveMonX$Sunlight, method = "spearman")##
## Spearman's rank correlation rho
##
## data: AveMonX$Precipitation and AveMonX$Sunlight
## S = 453, p-value = 0.047
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.58246cor.test(AveLomX$Precipitation, AveLomX$Irradiation, method = "pearson")##
## Pearson's product-moment correlation
##
## data: AveLomX$Precipitation and AveLomX$Irradiation
## t = -2.55, df = 10, p-value = 0.029
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.883188 -0.083241
## sample estimates:
## cor
## -0.62718cor.test(AveLomX$Precipitation, AveLomX$Irradiation, method = "spearman")##
## Spearman's rank correlation rho
##
## data: AveLomX$Precipitation and AveLomX$Irradiation
## S = 452, p-value = 0.052
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.58042cor.test(AveMonX$Precipitation, AveMonX$Irradiation, method = "pearson")##
## Pearson's product-moment correlation
##
## data: AveMonX$Precipitation and AveMonX$Irradiation
## t = -2.33, df = 10, p-value = 0.042
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.870921 -0.030231
## sample estimates:
## cor
## -0.59383cor.test(AveMonX$Precipitation, AveMonX$Irradiation, method = "spearman")##
## Spearman's rank correlation rho
##
## data: AveMonX$Precipitation and AveMonX$Irradiation
## S = 454, p-value = 0.049
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.58741Marvelous. Now, I need this to be separaed by habitat so I can
illustrate it. First, we will use the select function from
the dplyr package to specify the columns we wish to keep in
our new dataframes, Lome4 and Moun4. We will
then use the aggregate function to calculate means by Month
(Sample) and Age Group (Succession), creating
Lome4Ave and Moun4Ave. After renaming with
colnames, we merge with SolarMeltAve so that
we have values of sunlight in our new dataframes as well. We do this
once more with the Irad dataframe to include irradiation.
Note that this only works because all these dataframes are ordered in
the same way.
Lome4 <- Lomerio %>% dplyr::select(Sample, LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Detritus.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Total.Mg.ha.month, Precipitation, Succession)
Moun4 <- Mountain %>% dplyr::select(Sample, LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Detritus.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Total.Mg.ha.month, Precipitation, Succession)
Lome4Ave <- aggregate(Lome4[ , 2:8], by = list(Lome4$Sample, Lome4$Succession), FUN = mean)
Moun4Ave <- aggregate(Moun4[ , 2:8], by = list(Moun4$Sample, Moun4$Succession), FUN = mean, na.rm = TRUE)
colnames(Lome4Ave) <- c("Month", "Habitat", "LeafLitter.Mg.ha.month",
"Branches.Mg.ha.month", "Detritus.Mg.ha.month",
"Flowers.Mg.ha.month", "Fruits.Mg.ha.month",
"Total.Mg.ha.month", "Precipitation")
colnames(Moun4Ave) <- c("Month", "Habitat", "LeafLitter.Mg.ha.month",
"Branches.Mg.ha.month", "Detritus.Mg.ha.month",
"Flowers.Mg.ha.month", "Fruits.Mg.ha.month",
"Total.Mg.ha.month", "Precipitation")
summaryAveLom <- merge(SolarMeltAve, Lome4Ave, by="Month")
summaryAveMon <- merge(SolarMeltAve, Moun4Ave, by="Month")
summaryAveLomX <- merge(summaryAveLom, Irad, by="Month")
summaryAveMonX <- merge(summaryAveMon, Irad, by="Month") Next step is to set the Months into the proper format for figure
making. Once relabeled, we can use the as.Date function to
turn our variable into Date class. Then, we use the month
function to make it so that our date class variable shows only the
month. Finally, we rename the Age Groups into their final names.
levels(summaryAveLomX$Month) <- list("2017-08-01" = "Aug", "2017-09-01" = "Sept", "2017-10-01" = "Oct",
"2017-11-01" = "Nov", "2017-12-01" = "Dec", "2018-01-01" = "Jan",
"2018-02-01" = "Feb", "2018-03-01" = "Mar", "2018-04-01" = "Apr",
"2018-05-01" = "May", "2018-06-01" = "Jun","2018-07-01" = "Jul")
levels(summaryAveMonX$Month) <- list("2017-09-01" = "Sept", "2017-10-01" = "Oct",
"2017-11-01" = "Nov", "2017-12-01" = "Dec", "2018-01-01" = "Jan",
"2018-02-01" = "Feb", "2018-03-01" = "Mar", "2018-04-01" = "Apr",
"2018-05-01" = "May", "2018-06-01" = "Jun","2018-07-01" = "Jul",
"2018-08-01" = "Aug")
summaryAveLomX$Month <- as.Date(summaryAveLomX$Month)
summaryAveLomX$Month <- month(summaryAveLomX$Month, label = TRUE)
summaryAveMonX$Month <- as.Date(summaryAveMonX$Month)
summaryAveMonX$Month <- month(summaryAveMonX$Month, label = TRUE)
levels(summaryAveLomX$Habitat) <- list(EH = "DL1", IH = "DL2", OH = "RL")
levels(summaryAveMonX$Habitat) <- list(EM = "DM1", IM = "DM2", OM = "RM")We’re almost ready to graph! However, this complicated graph has three y-axes. I can only in R insert two, so we’ll have to make two versions of this graph that we at the end combine in PowerPoint.
Figure 4
This first version will include only the Precipitation as the secondary y-axis. We will create a graph using ggplot for each landscape and store them before combining them with a shared legend.
GLX1 <- ggplot(data=summaryAveLomX, aes(x=Month, y=LeafLitter.Mg.ha.month )) +
geom_bar(stat="identity", aes(fill = Habitat), position="dodge", alpha = 0.8) +
xlab("") +
ylab(bquote('Leaf Litter '(Mg~ha^-1~mth^-1))) +
geom_line(aes(Month, Irradiation*.003, group = Habitat), color = "GoldenRod") +
geom_point(aes(Month, Irradiation*.003, group = Habitat), color = "Black", shape = "☼", size = 3) +
geom_boxplot(aes(Month, Precipitation*.002), color = "darkblue") +
scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1),
sec.axis = sec_axis(~./.002, name="Precipitation (mm)", breaks = seq(0,500,100))) +
scale_fill_manual(values=c("#82C27B", "#43833C", "#23431F")) +
expand_limits(y = c(0, 1.1)) +
guides(fill=guide_legend(title="Age Group")) +
theme_bw(base_size = 12) + theme(legend.title = element_text(color = "black"))
GMX1 <- ggplot(data=summaryAveMonX, aes(x=Month, y=LeafLitter.Mg.ha.month )) +
geom_bar(stat="identity", aes(fill = Habitat), position="dodge", alpha = 0.8) +
xlab("Month") +
ylab(bquote('Leaf Litter '(Mg~ha^-1~mth^-1))) +
geom_line(aes(Month, Irradiation*.003, group = Habitat), color = "GoldenRod") +
geom_point(aes(Month, Irradiation*.003, group = Habitat), color = "Black", shape = "☼", size = 3) +
geom_boxplot(aes(Month, Precipitation*.002), color = "darkblue") +
scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1),
sec.axis = sec_axis(~./0.002, name="Precipitation (mm)", breaks = seq(0,500,100))) +
scale_fill_manual(values=c("#FBBA84", "#CE803F", "#67360D")) +
expand_limits(y = c(0, 1.1)) +
guides(fill=guide_legend(title="")) +
theme(legend.key.size = unit(0.2, 'cm')) +
theme_bw(base_size = 12)
combined1 <- GLX1 + GMX1 & theme(legend.position = "bottom") +
theme(legend.title = element_text(color = "black")) + theme(legend.title.align=0.5) +
theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'))
combined1 + plot_layout(guides = "collect", nrow = 2) + plot_annotation(tag_levels = c('a'))
GLX2 <- ggplot(data=summaryAveLomX, aes(x=Month, y=LeafLitter.Mg.ha.month )) +
geom_bar(stat="identity", aes(fill = Habitat), position="dodge", alpha = 0.8) +
xlab("") +
ylab(bquote('Leaf Litter '(Mg~ha^-1~mth^-1))) +
geom_line(aes(Month, Irradiation*0.003, group = Habitat), color = "GoldenRod") +
geom_point(aes(Month, Irradiation*0.003, group = Habitat), color = "Black", shape = "☼", size = 3) +
geom_boxplot(aes(Month, Precipitation*0.002), color = "darkblue") +
scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1),
sec.axis = sec_axis(~./0.003, name= bquote('Irradiation'~(W~m^-2)))) +
scale_fill_manual(values=c("#82C27B", "#43833C", "#23431F")) +
expand_limits(y = c(0, 1.1)) +
guides(fill=guide_legend(title="Age Group")) +
theme_bw(base_size = 12) + theme(legend.title = element_text(color = "black"))
GMX2 <- ggplot(data=summaryAveMonX, aes(x=Month, y=LeafLitter.Mg.ha.month )) +
geom_bar(stat="identity", aes(fill = Habitat), position="dodge", alpha = 0.8) +
xlab("Month") +
ylab(bquote('Leaf Litter '(Mg~ha^-1~mth^-1))) +
geom_line(aes(Month, Irradiation*0.003, group = Habitat), color = "GoldenRod") +
geom_point(aes(Month, Irradiation*0.003, group = Habitat), color = "Black", shape = "☼", size = 3) +
geom_boxplot(aes(Month, Precipitation*0.002), color = "darkblue") +
scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1),
sec.axis = sec_axis(~./0.003, name= bquote('Irradiation'~(W~m^-2)))) +
scale_fill_manual(values=c("#FBBA84", "#CE803F", "#67360D")) +
expand_limits(y = c(0, 1.1)) +
guides(fill=guide_legend(title="")) +
theme_bw(base_size = 12) + theme(legend.title = element_text(color = "black"))
combined2 <- GLX2 + GMX2 & theme(legend.position = "bottom") +
theme(legend.title = element_text(color = "black")) + theme(legend.title.align=0.5) +
theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'))
combined2 + plot_layout(guides = "collect", nrow = 2) + plot_annotation(tag_levels = c('a'))
Statistical Analysis
Litter x Landscape
Does litterfall differ between lomerio and mountain? We can find out through Analysis of Variance (ANOVA) available through the anova function. We will use the appropriate transformations for each variable. For interpretation, typically a p-value below 0.05 is considered statistically significant, and a p-value between 0.05 and 0.1 as marginally significant.
anova(lm(data = LeafLitt2, sqrt(LeafLitter.Mg.ha.month) ~ Paisaje,na.action=na.exclude))## Analysis of Variance Table
##
## Response: sqrt(LeafLitter.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Paisaje 1 1.29 1.292 47 3e-11 ***
## Residuals 367 10.09 0.027
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(data = LeafLitt2, Branches.Mg.ha.month^1/3 ~ Paisaje,na.action=na.exclude))## Analysis of Variance Table
##
## Response: Branches.Mg.ha.month^1/3
## Df Sum Sq Mean Sq F value Pr(>F)
## Paisaje 1 0.0012 0.001178 7.51 0.0064 **
## Residuals 367 0.0575 0.000157
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(data = LeafLitt2, log(Fruits.Mg.ha.month+0.05) ~ Paisaje,na.action=na.exclude))## Analysis of Variance Table
##
## Response: log(Fruits.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## Paisaje 1 0.2 0.156 1.7 0.19
## Residuals 367 33.7 0.092anova(lm(data = LeafLitt2, log(Flowers.Mg.ha.month+0.05) ~ Paisaje,na.action=na.exclude))## Analysis of Variance Table
##
## Response: log(Flowers.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## Paisaje 1 0.05 0.0509 4.48 0.035 *
## Residuals 367 4.17 0.0114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(data = LeafLitt2, sqrt(Detritus.Mg.ha.month) ~ Paisaje,na.action=na.exclude))## Analysis of Variance Table
##
## Response: sqrt(Detritus.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Paisaje 1 0.00 0.00017 0.01 0.91
## Residuals 367 4.82 0.01314anova(lm(data = LeafLitt2, sqrt(Total.Mg.ha.month) ~ Paisaje,na.action=na.exclude))## Analysis of Variance Table
##
## Response: sqrt(Total.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Paisaje 1 0.84 0.842 23.3 2.1e-06 ***
## Residuals 367 13.28 0.036
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Figure 2
Lets make a stacked barplot to visualize how leaf litter is composed
in each landscape. To do this, first we need to make a new dataframe,
LomerioStack1, that will hold only the columsn that we are
interested in. For this, we can use the -which function to
select the columns we DONT want by name. Yeah, we could have just used
the which function and select the ones we want, but, we’re
still learning!
LomerioStack1 <- Lomerio[ , -which(names(Lomerio) %in% c("Record","Paisaje", "Locality",
"Successional.State", "Plot",
"LeafLitter.g", "LeafLitter.Mg.ha.month2",
"Branches.g", "Detritus.g",
"Flowers.g", "Fruits.g", "Total.g",
"Total.Mg.ha.month", "Precipitation","Temp.Ave","Temp.Max",
"PlotSample", "CodeX", "CodeZ","CodeY", "RainNow",
"Rain1", "Rain2", "Rain3", "TempNow",
"Temp1", "Temp2", "Temp3"))]Next we use the melt function from the
reshape package to change the format of our table so that
we have a new column that holds all litter types. We’ll call this new
dataframe LomerioStack2. We will then use
colnames to rename the columns into what we want. Next, we
use the summarise function from the dplyr
package to calculate mean values of litter production for each specific
grouping. Why mean and not sum? Well, because our landscapes do not have
an equal number of plots, and we’re more interested in the proportion of
littertype. We’ll throw those mean values into a new dataframe,
LomerioStack3.
LomerioStack2 <- melt(LomerioStack1, id.vars=c("Sample", "Succession"))
colnames(LomerioStack2) <- c("Month", "Habitat", "LitterType", "Value")
LomerioStack3 <- as.data.frame(LomerioStack2 %>%
group_by(LitterType, Habitat) %>%
dplyr::summarise(Value = mean(Value)))We’ll do the exact same for the mountain plots, with the exception
that halfway through we’ll use the na.omit function to
remove a few NA values.
MountainStack1 <- Mountain[ , -which(names(Mountain) %in% c("Record","Paisaje", "Locality", "Successional.State", "Plot", "LeafLitter.g", "LeafLitter.Mg.ha.month2",
"Branches.g", "Detritus.g", "Flowers.g", "Fruits.g", "Total.g",
"Total.Mg.ha.month", "Precipitation","Temp.Ave","Temp.Max", "PlotSample",
"CodeX", "CodeZ","CodeY", "RainNow", "Rain1", "Rain2", "Rain3", "TempNow",
"Temp1", "Temp2", "Temp3"))]
MountainStack2 <- melt(MountainStack1, id.vars=c("Sample", "Succession"))
colnames(MountainStack2) <- c("Month", "Habitat", "LitterType", "Value")
MountainStack2.5 <- na.omit(MountainStack2)
MountainStack3 <- as.data.frame(MountainStack2.5 %>%
group_by(LitterType, Habitat) %>%
dplyr::summarise(Value = mean(Value)))Time to graph! We’ll use ggplot to create this stacked
barplots, but first, we will rename the landscapes to their english
names and bind them together by row using the rbind
function, creating the dataframe Megastack. We’ll also
rename the Age Groups into the final names the publication needs.
Finally, with colnames, we rename the columns.
LomerioStack3$Landscape <- "Hill"
MountainStack3$Landscape <- "Mountain"
Megastack <- rbind(LomerioStack3, MountainStack3)
levels(Megastack$Habitat) <- list(EH = "DL1", IH = "DL2", OH = "RL",
EM = "DM1", IM = "DM2", OM = "RM")
colnames(Megastack) <- c("LitterType", "Age Group", "Value", "Landscape")ggplot(Megastack, aes(fill=LitterType, y=Value, x=`Age Group`)) +
geom_bar(position="stack", stat="identity", color = "black", alpha = 0.7) +
ylab(bquote('Mean Total Litter '(Mg~ha^-1~mth^-1))) +
ylim(0, 1) +
labs(fill = "") +
scale_fill_brewer(palette = "Accent", labels=c("Leaf", "Branch", "Detritus",
"Flower", "Fruit"),
guide = guide_legend(direction = "horizontal",
title.position = "left")) +
theme_bw(base_size = 13) + theme(legend.title = element_text(color = "black")) +
theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid')) +
theme(legend.position = "bottom") +
theme(legend.title.align=0.5) +
theme(legend.key.size = unit(0.4, 'cm')) +
facet_wrap(~Landscape, scales = "free") 
Litter x Precipitation
Do the habitats differ in their relationship between leaf litter and
precipitation? We can find out by running an ANCOVA with
Precipitation and Habitat as interactive terms
in an anova function.
anova(lm(data = summaryAveLom, sqrt(LeafLitter.Mg.ha.month) ~ Precipitation * Habitat))## Analysis of Variance Table
##
## Response: sqrt(LeafLitter.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Precipitation 1 0.0344 0.0344 3.94 0.056 .
## Habitat 2 0.2788 0.1394 15.98 1.9e-05 ***
## Precipitation:Habitat 2 0.0030 0.0015 0.17 0.841
## Residuals 30 0.2616 0.0087
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(data = summaryAveMon, sqrt(LeafLitter.Mg.ha.month) ~ Precipitation * Habitat))## Analysis of Variance Table
##
## Response: sqrt(LeafLitter.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Precipitation 1 0.1527 0.1527 24.23 2.9e-05 ***
## Habitat 2 0.0789 0.0395 6.26 0.0053 **
## Precipitation:Habitat 2 0.0025 0.0012 0.20 0.8238
## Residuals 30 0.1890 0.0063
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1They do not, so we will be pooling the habitats within landscape together.
Here we’re going to include interactions between Precipitation and Month as well as Habitat and Month.
Hill
anova(lm(sqrt(LeafLitter.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: sqrt(LeafLitter.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.923 0.461 25.43 6.8e-10 ***
## RainNow 1 0.126 0.126 6.94 0.0096 **
## Sample 11 0.936 0.085 4.69 6.4e-06 ***
## Succession:RainNow 2 0.008 0.004 0.21 0.8076
## Succession:Sample 22 0.351 0.016 0.88 0.6221
## Residuals 117 2.122 0.018
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(LeafLitter.Mg.ha.month) ~ Succession + RainNow + Sample + Succession + Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: sqrt(LeafLitter.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.923 0.461 26.22 2.1e-10 ***
## RainNow 1 0.126 0.126 7.15 0.0084 **
## Sample 11 0.936 0.085 4.84 2.6e-06 ***
## Residuals 141 2.481 0.018
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(Branches.Mg.ha.month^(1/3) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: Branches.Mg.ha.month^(1/3)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.220 0.1102 16.70 4.2e-07 ***
## RainNow 1 0.013 0.0131 1.99 0.16126
## Sample 11 0.270 0.0245 3.72 0.00014 ***
## Succession:RainNow 2 0.002 0.0012 0.17 0.83982
## Succession:Sample 22 0.124 0.0056 0.85 0.65681
## Residuals 117 0.772 0.0066
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(Branches.Mg.ha.month^(1/3) ~ Succession + RainNow + Sample + Succession + Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: Branches.Mg.ha.month^(1/3)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.220 0.1102 17.30 1.9e-07 ***
## RainNow 1 0.013 0.0131 2.06 0.15
## Sample 11 0.270 0.0245 3.86 7.2e-05 ***
## Residuals 141 0.898 0.0064
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(log(Flowers.Mg.ha.month+0.05) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: log(Flowers.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.068 0.0338 5.61 0.0047 **
## RainNow 1 0.001 0.0005 0.08 0.7724
## Sample 11 0.134 0.0122 2.02 0.0320 *
## Succession:RainNow 2 0.028 0.0138 2.29 0.1060
## Succession:Sample 22 0.247 0.0112 1.86 0.0185 *
## Residuals 117 0.706 0.0060
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(log(Flowers.Mg.ha.month+0.05) ~ Succession + RainNow + Sample + Succession + Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: log(Flowers.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.068 0.0338 4.87 0.0091 **
## RainNow 1 0.001 0.0005 0.07 0.7876
## Sample 11 0.134 0.0122 1.76 0.0675 .
## Residuals 141 0.980 0.0070
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(log(Fruits.Mg.ha.month+0.05) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: log(Fruits.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.33 0.1662 1.61 0.21
## RainNow 1 0.09 0.0869 0.84 0.36
## Sample 11 1.94 0.1765 1.71 0.08 .
## Succession:RainNow 2 0.23 0.1171 1.13 0.33
## Succession:Sample 22 1.64 0.0746 0.72 0.81
## Residuals 117 12.11 0.1035
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(log(Fruits.Mg.ha.month+0.05) ~ Succession+ RainNow + Sample + Succession + Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: log(Fruits.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.33 0.1662 1.68 0.191
## RainNow 1 0.09 0.0869 0.88 0.351
## Sample 11 1.94 0.1765 1.78 0.063 .
## Residuals 141 13.98 0.0992
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Detritus.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: sqrt(Detritus.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.311 0.1555 17.57 2.1e-07 ***
## RainNow 1 0.104 0.1037 11.71 0.00086 ***
## Sample 11 0.131 0.0119 1.34 0.20902
## Succession:RainNow 2 0.052 0.0260 2.94 0.05676 .
## Succession:Sample 22 0.196 0.0089 1.01 0.46200
## Residuals 117 1.035 0.0089
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Detritus has marginally sig interaction Succession:RainNow.
anova(lm(sqrt(Total.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: sqrt(Total.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 1.382 0.691 31.71 9.9e-12 ***
## RainNow 1 0.236 0.236 10.84 0.00131 **
## Sample 11 0.892 0.081 3.72 0.00014 ***
## Succession:RainNow 2 0.027 0.013 0.61 0.54305
## Succession:Sample 22 0.451 0.020 0.94 0.54392
## Residuals 117 2.549 0.022
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Total.Mg.ha.month) ~ Succession+ RainNow + Sample + Succession + Sample, data = Lomerio))## Analysis of Variance Table
##
## Response: sqrt(Total.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 1.382 0.691 32.19 3.1e-12 ***
## RainNow 1 0.236 0.236 11.00 0.0012 **
## Sample 11 0.892 0.081 3.78 9.3e-05 ***
## Residuals 141 3.027 0.021
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Because only Detritus had a marginally significant interaction, for
detritus only we will separate by Age Group and evaluate the
relationship between rain and detritus for each group. First, we use
droplevels and subset to create these groups.
Then, we run the anovas again.
LomeDL1 <- droplevels(subset(Lomerio, Lomerio$Succession == "DL1"))
LomeDL2 <- droplevels(subset(Lomerio, Lomerio$Succession == "DL2"))
LomeRL <- droplevels(subset(Lomerio, Lomerio$Succession == "RL"))
anova(lm(sqrt(Detritus.Mg.ha.month) ~ RainNow + Sample, data = LomeDL1)) ## Analysis of Variance Table
##
## Response: sqrt(Detritus.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## RainNow 1 0.0033 0.00334 0.51 0.48
## Sample 11 0.0847 0.00770 1.17 0.36
## Residuals 23 0.1517 0.00659anova(lm(sqrt(Detritus.Mg.ha.month) ~ RainNow + Sample, data = LomeDL2)) ## Analysis of Variance Table
##
## Response: sqrt(Detritus.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## RainNow 1 0.097 0.0973 10.33 0.0021 **
## Sample 11 0.129 0.0117 1.25 0.2775
## Residuals 59 0.555 0.0094
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Detritus.Mg.ha.month) ~ RainNow + Sample, data = LomeRL)) ## Analysis of Variance Table
##
## Response: sqrt(Detritus.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## RainNow 1 0.057 0.0569 6.06 0.019 *
## Sample 11 0.111 0.0101 1.08 0.407
## Residuals 35 0.329 0.0094
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Next we use the cor.test function to obtain Pearson’s and Spearman’s Rank coefficients.
cor.test(sqrt(Lomerio$LeafLitter.Mg.ha.month), Lomerio$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: sqrt(Lomerio$LeafLitter.Mg.ha.month) and Lomerio$RainNow
## t = -2.07, df = 154, p-value = 0.041
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3132760 -0.0072663
## sample estimates:
## cor
## -0.16422cor.test(Lomerio$LeafLitter.Mg.ha.month, Lomerio$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Lomerio$LeafLitter.Mg.ha.month and Lomerio$RainNow
## S = 727198, p-value = 0.063
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.14934cor.test(Lomerio$Branches.Mg.ha.month^(1/3), Lomerio$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: Lomerio$Branches.Mg.ha.month^(1/3) and Lomerio$RainNow
## t = -1.16, df = 154, p-value = 0.25
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.246880 0.064729
## sample estimates:
## cor
## -0.093361cor.test(Lomerio$Branches.Mg.ha.month, Lomerio$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Lomerio$Branches.Mg.ha.month and Lomerio$RainNow
## S = 703803, p-value = 0.16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.11236cor.test(log(Lomerio$Flowers.Mg.ha.month+0.05), Lomerio$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: log(Lomerio$Flowers.Mg.ha.month + 0.05) and Lomerio$RainNow
## t = -0.237, df = 154, p-value = 0.81
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.17569 0.13848
## sample estimates:
## cor
## -0.019076cor.test(Lomerio$Flowers.Mg.ha.month, Lomerio$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Lomerio$Flowers.Mg.ha.month and Lomerio$RainNow
## S = 660023, p-value = 0.59
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.043168cor.test(log(Lomerio$Fruits.Mg.ha.month+0.05), Lomerio$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: log(Lomerio$Fruits.Mg.ha.month + 0.05) and Lomerio$RainNow
## t = -0.892, df = 154, p-value = 0.37
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.226265 0.086442
## sample estimates:
## cor
## -0.071672cor.test(Lomerio$Fruits.Mg.ha.month, Lomerio$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Lomerio$Fruits.Mg.ha.month and Lomerio$RainNow
## S = 682877, p-value = 0.33
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.079289cor.test(sqrt(Lomerio$Total.Mg.ha.month), Lomerio$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: sqrt(Lomerio$Total.Mg.ha.month) and Lomerio$RainNow
## t = -2.57, df = 154, p-value = 0.011
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.348510 -0.046806
## sample estimates:
## cor
## -0.20246cor.test(Lomerio$Total.Mg.ha.month, Lomerio$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Lomerio$Total.Mg.ha.month and Lomerio$RainNow
## S = 742872, p-value = 0.03
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.17411cor.test(sqrt(LomeDL1$Detritus.Mg.ha.month), LomeDL1$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: sqrt(LomeDL1$Detritus.Mg.ha.month) and LomeDL1$RainNow
## t = 0.693, df = 34, p-value = 0.49
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.21905 0.42985
## sample estimates:
## cor
## 0.11797cor.test(LomeDL1$Detritus.Mg.ha.month, LomeDL1$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: LomeDL1$Detritus.Mg.ha.month and LomeDL1$RainNow
## S = 7450, p-value = 0.81
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.041189cor.test(sqrt(LomeDL2$Detritus.Mg.ha.month), LomeDL2$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: sqrt(LomeDL2$Detritus.Mg.ha.month) and LomeDL2$RainNow
## t = -3.15, df = 70, p-value = 0.0024
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.54026 -0.13185
## sample estimates:
## cor
## -0.35274cor.test(LomeDL2$Detritus.Mg.ha.month, LomeDL2$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: LomeDL2$Detritus.Mg.ha.month and LomeDL2$RainNow
## S = 82770, p-value = 0.0045
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.3308cor.test(sqrt(LomeRL$Detritus.Mg.ha.month), LomeRL$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: sqrt(LomeRL$Detritus.Mg.ha.month) and LomeRL$RainNow
## t = -2.44, df = 46, p-value = 0.019
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.56804 -0.06021
## sample estimates:
## cor
## -0.33855cor.test(LomeRL$Detritus.Mg.ha.month, LomeRL$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: LomeRL$Detritus.Mg.ha.month and LomeRL$RainNow
## S = 25976, p-value = 0.0038
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.40988Mountain
We’ll repeat the process for the mountain landscape.
anova(lm(sqrt(LeafLitter.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: sqrt(LeafLitter.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.531 0.266 15.03 9.5e-07 ***
## RainNow 1 0.781 0.781 44.17 3.7e-10 ***
## Sample 11 0.881 0.080 4.53 5.0e-06 ***
## Succession:RainNow 2 0.006 0.003 0.17 0.84
## Succession:Sample 22 0.349 0.016 0.90 0.60
## Residuals 174 3.075 0.018
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(LeafLitter.Mg.ha.month) ~ Succession + RainNow + Sample + Succession + Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: sqrt(LeafLitter.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.53 0.266 15.34 6.4e-07 ***
## RainNow 1 0.78 0.781 45.07 2.0e-10 ***
## Sample 11 0.88 0.080 4.62 2.9e-06 ***
## Residuals 198 3.43 0.017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(Branches.Mg.ha.month^(1/3) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: Branches.Mg.ha.month^(1/3)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.328 0.1638 21.34 5.1e-09 ***
## RainNow 1 0.001 0.0015 0.19 0.66
## Sample 11 0.313 0.0285 3.71 9.2e-05 ***
## Succession:RainNow 2 0.006 0.0031 0.40 0.67
## Succession:Sample 22 0.146 0.0066 0.86 0.64
## Residuals 174 1.335 0.0077
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(Branches.Mg.ha.month^(1/3) ~ Succession+ RainNow + Sample + Succession + Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: Branches.Mg.ha.month^(1/3)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.328 0.1638 21.80 2.8e-09 ***
## RainNow 1 0.001 0.0015 0.20 0.66
## Sample 11 0.313 0.0285 3.79 6.0e-05 ***
## Residuals 198 1.487 0.0075
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(log(Flowers.Mg.ha.month+0.05) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: log(Flowers.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.119 0.0593 4.94 0.0082 **
## RainNow 1 0.001 0.0009 0.07 0.7849
## Sample 11 0.218 0.0198 1.65 0.0897 .
## Succession:RainNow 2 0.060 0.0298 2.48 0.0871 .
## Succession:Sample 22 0.501 0.0228 1.89 0.0125 *
## Residuals 174 2.091 0.0120
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(log(Fruits.Mg.ha.month+0.05) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: log(Fruits.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 1.80 0.902 13.22 4.5e-06 ***
## RainNow 1 0.20 0.202 2.96 0.0871 .
## Sample 11 2.26 0.205 3.01 0.0011 **
## Succession:RainNow 2 0.21 0.104 1.53 0.2195
## Succession:Sample 22 1.06 0.048 0.70 0.8313
## Residuals 174 11.87 0.068
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(log(Fruits.Mg.ha.month+0.05) ~ Succession+ RainNow + Sample + Succession + Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: log(Fruits.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 1.80 0.902 13.60 2.9e-06 ***
## RainNow 1 0.20 0.202 3.04 0.08257 .
## Sample 11 2.26 0.205 3.09 0.00074 ***
## Residuals 198 13.14 0.066
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Detritus.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: sqrt(Detritus.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.533 0.2663 25.48 2e-10 ***
## RainNow 1 0.013 0.0129 1.24 0.26728
## Sample 11 0.362 0.0329 3.15 0.00066 ***
## Succession:RainNow 2 0.010 0.0052 0.50 0.60784
## Succession:Sample 22 0.255 0.0116 1.11 0.34213
## Residuals 174 1.819 0.0105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Detritus.Mg.ha.month) ~ Succession+ RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: sqrt(Detritus.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.533 0.2663 25.15 2.5e-10 ***
## RainNow 1 0.013 0.0129 1.22 0.27033
## Sample 11 0.362 0.0329 3.11 0.00076 ***
## Succession:Sample 22 0.221 0.0100 0.95 0.53420
## Residuals 176 1.863 0.0106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Total.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: sqrt(Total.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 1.33 0.665 29.43 9.8e-12 ***
## RainNow 1 0.59 0.586 25.94 9.1e-07 ***
## Sample 11 1.33 0.121 5.34 2.8e-07 ***
## Succession:RainNow 2 0.00 0.001 0.06 0.94
## Succession:Sample 22 0.56 0.025 1.12 0.33
## Residuals 174 3.93 0.023
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Total.Mg.ha.month) ~ Succession+ RainNow + Sample + Succession + Sample, data = Mountain, na.action=na.exclude))## Analysis of Variance Table
##
## Response: sqrt(Total.Mg.ha.month)
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 1.33 0.665 29.32 7.1e-12 ***
## RainNow 1 0.59 0.586 25.84 8.6e-07 ***
## Sample 11 1.33 0.121 5.32 2.3e-07 ***
## Residuals 198 4.49 0.023
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.test(sqrt(Mountain$LeafLitter.Mg.ha.month), Mountain$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: sqrt(Mountain$LeafLitter.Mg.ha.month) and Mountain$RainNow
## t = -5.76, df = 211, p-value = 3e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.47906 -0.24611
## sample estimates:
## cor
## -0.36835cor.test(Mountain$LeafLitter.Mg.ha.month, Mountain$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Mountain$LeafLitter.Mg.ha.month and Mountain$RainNow
## S = 2164663, p-value = 2.6e-07
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.34404cor.test(Mountain$Branches.Mg.ha.month^(1/3), Mountain$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: Mountain$Branches.Mg.ha.month^(1/3) and Mountain$RainNow
## t = -0.335, df = 211, p-value = 0.74
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.15700 0.11172
## sample estimates:
## cor
## -0.023054cor.test(Mountain$Branches.Mg.ha.month, Mountain$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Mountain$Branches.Mg.ha.month and Mountain$RainNow
## S = 1600678, p-value = 0.93
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.0061382cor.test(log(Mountain$Flowers.Mg.ha.month+0.05), Mountain$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: log(Mountain$Flowers.Mg.ha.month + 0.05) and Mountain$RainNow
## t = 0.315, df = 211, p-value = 0.75
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.11311 0.15563
## sample estimates:
## cor
## 0.021655cor.test(Mountain$Flowers.Mg.ha.month, Mountain$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Mountain$Flowers.Mg.ha.month and Mountain$RainNow
## S = 1674001, p-value = 0.57
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.039388cor.test(log(Mountain$Fruits.Mg.ha.month+0.05), Mountain$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: log(Mountain$Fruits.Mg.ha.month + 0.05) and Mountain$RainNow
## t = 1.65, df = 211, p-value = 0.1
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.021708 0.243780
## sample estimates:
## cor
## 0.11305cor.test(Mountain$Fruits.Mg.ha.month, Mountain$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Mountain$Fruits.Mg.ha.month and Mountain$RainNow
## S = 1456977, p-value = 0.17
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.095362cor.test(sqrt(Mountain$Detritus.Mg.ha.month), Mountain$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: sqrt(Mountain$Detritus.Mg.ha.month) and Mountain$RainNow
## t = -0.908, df = 211, p-value = 0.37
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.195169 0.072667
## sample estimates:
## cor
## -0.062374cor.test(Mountain$Detritus.Mg.ha.month, Mountain$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Mountain$Detritus.Mg.ha.month and Mountain$RainNow
## S = 1700554, p-value = 0.42
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.055875cor.test(sqrt(Mountain$Total.Mg.ha.month), Mountain$RainNow, method = "pearson")##
## Pearson's product-moment correlation
##
## data: sqrt(Mountain$Total.Mg.ha.month) and Mountain$RainNow
## t = -4.08, df = 211, p-value = 6.3e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.39077 -0.14126
## sample estimates:
## cor
## -0.27055cor.test(Mountain$Total.Mg.ha.month, Mountain$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: Mountain$Total.Mg.ha.month and Mountain$RainNow
## S = 2018308, p-value = 0.00019
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.25317We had a marginally significant interaction with flowers, so lets separate.
MounDM1 <- droplevels(subset(Mountain, Mountain$Succession == "DM1"))
MounDM2 <- droplevels(subset(Mountain, Mountain$Succession == "DM2"))
MounRM <- droplevels(subset(Mountain, Mountain$Succession == "RM"))
anova(lm(log(Flowers.Mg.ha.month+0.05) ~ RainNow + Sample, data = MounDM1)) # 0.08966## Analysis of Variance Table
##
## Response: log(Flowers.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## RainNow 1 0.016 0.01597 2.98 0.09 .
## Sample 11 0.089 0.00808 1.51 0.15
## Residuals 59 0.317 0.00536
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(log(Flowers.Mg.ha.month+0.05) ~ RainNow + Sample, data = MounDM2)) # 0.06538## Analysis of Variance Table
##
## Response: log(Flowers.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## RainNow 1 0.034 0.0342 3.53 0.065 .
## Sample 11 0.093 0.0085 0.88 0.568
## Residuals 58 0.562 0.0097
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(log(Flowers.Mg.ha.month+0.05) ~ RainNow + Sample, data = MounRM)) # 0.98921## Analysis of Variance Table
##
## Response: log(Flowers.Mg.ha.month + 0.05)
## Df Sum Sq Mean Sq F value Pr(>F)
## RainNow 1 0.000 0.0000 0.00 0.989
## Sample 11 0.546 0.0497 2.34 0.019 *
## Residuals 57 1.212 0.0213
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.test(MounDM1$Flowers.Mg.ha.month, MounDM1$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: MounDM1$Flowers.Mg.ha.month and MounDM1$RainNow
## S = 72221, p-value = 0.18
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.16119cor.test(MounDM2$Flowers.Mg.ha.month, MounDM2$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: MounDM2$Flowers.Mg.ha.month and MounDM2$RainNow
## S = 59900, p-value = 0.97
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.0043576cor.test(MounRM$Flowers.Mg.ha.month, MounRM$RainNow, method = "spearman")##
## Spearman's rank correlation rho
##
## data: MounRM$Flowers.Mg.ha.month and MounRM$RainNow
## S = 56523, p-value = 0.93
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.011066Litter x Structure
Do the litter traits correlate with the tree community variables?
Before we’re able to answer this question, we need to calculate total
litter and precipitation by our CodeZ variable. We can
accomplish this with the aggregate function. We’ll then
merge our dataframes together and rename using
colnames.
DFX1 <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month,
Flowers.Mg.ha.month, Fruits.Mg.ha.month,
Detritus.Mg.ha.month, Total.Mg.ha.month) ~ CodeZ,
data = LeafLitt2, FUN = sum, na.action=na.exclude)
DFX2 <- aggregate(cbind(RainNow, Rain1, Rain2, Rain3) ~ CodeZ,
data = LeafLitt2, FUN = mean, na.action=na.exclude)
DFX3 <- merge(DFX1, DFX2, by="CodeZ")
colnames(DFX3) <- c("CodeZ", "Leaf","Branch","Flower","Fruit","Detritus","Total",
"RainNow","Rain1","Rain2","Rain3")Next we merge our new dataframe, DFX3, with the
Struc dataframe that contains basal area, stem density,
species density, biomass, and species diversity. We then
subset to create a dataframe for Hill and one for Mountain.
We’ll also subset again to create dataframes for specific age groups
within landscape.
Struc <- merge(Struc, DFX3, by="CodeZ") # After merging we split by landscape.
StrucLome <- droplevels(subset(Struc, Struc$Landscape == "Lomerio"))
StrucMon <- droplevels(subset(Struc, Struc$Landscape == "Mountain"))
StrucDL1 <- droplevels(subset(Struc, Struc$Succession == "DL1")) # Lomerio
StrucDL2 <- droplevels(subset(Struc, Struc$Succession == "DL2"))
StrucRL <- droplevels(subset(Struc, Struc$Succession == "RL"))
StrucDM1 <- droplevels(subset(Struc, Struc$Succession == "DM1")) # Mountain
StrucDM2 <- droplevels(subset(Struc, Struc$Succession == "DM2"))
StrucRM <- droplevels(subset(Struc, Struc$Succession == "RM"))Now we can look at correlations. We will use the non-parametric Spearman’s Rank Correlation for this.
Leaf Litter
First, we’ll check out the correlations without separating by Age Group. We’ll start with leaf litter.
cor.test(Struc$Leaf, Struc$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Leaf and Struc$TreesHA
## S = 4735, p-value = 0.81
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.045372cor.test(Struc$Leaf, Struc$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Leaf and Struc$BasalHA
## S = 3743, p-value = 0.18
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.24531cor.test(Struc$Leaf, Struc$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Leaf and Struc$Biomass
## S = 3718, p-value = 0.17
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.25043cor.test(Struc$Leaf, Struc$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Leaf and Struc$Richness
## S = 3615, p-value = 0.14
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.27127cor.test(Struc$Leaf, Struc$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Leaf and Struc$Simpson_1.D
## S = 3362, p-value = 0.077
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.32226Next we will use ANCOVAs to determine if landscapes differ in their relationship between leaf litter and the structural variables.
anova(lm(sqrt(Leaf) ~ TreesHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Leaf)
## Df Sum Sq Mean Sq F value Pr(>F)
## TreesHA 1 0.083 0.083 1.68 0.2062
## Landscape 1 1.730 1.730 35.04 2.6e-06 ***
## TreesHA:Landscape 1 0.590 0.590 11.95 0.0018 **
## Residuals 27 1.333 0.049
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Leaf) ~ BasalHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Leaf)
## Df Sum Sq Mean Sq F value Pr(>F)
## BasalHA 1 0.226 0.226 4.01 0.0554 .
## Landscape 1 1.501 1.501 26.59 2e-05 ***
## BasalHA:Landscape 1 0.486 0.486 8.61 0.0068 **
## Residuals 27 1.524 0.056
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Leaf) ~ Biomass * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Leaf)
## Df Sum Sq Mean Sq F value Pr(>F)
## Biomass 1 0.269 0.269 4.59 0.041 *
## Landscape 1 1.460 1.460 24.97 3.1e-05 ***
## Biomass:Landscape 1 0.429 0.429 7.35 0.012 *
## Residuals 27 1.578 0.058
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Leaf) ~ Richness * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Leaf)
## Df Sum Sq Mean Sq F value Pr(>F)
## Richness 1 0.315 0.315 5.16 0.031 *
## Landscape 1 1.412 1.412 23.16 5e-05 ***
## Richness:Landscape 1 0.362 0.362 5.94 0.022 *
## Residuals 27 1.647 0.061
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Leaf) ~ Simpson_1.D * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Leaf)
## Df Sum Sq Mean Sq F value Pr(>F)
## Simpson_1.D 1 0.231 0.231 2.77 0.10735
## Landscape 1 1.234 1.234 14.82 0.00066 ***
## Simpson_1.D:Landscape 1 0.021 0.021 0.25 0.62019
## Residuals 27 2.250 0.083
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There was a significant interaction in all but diversity. Next we separate by landscape.
cor.test(StrucLome$Leaf, StrucLome$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Leaf and StrucLome$TreesHA
## S = 62, p-value = 0.00078
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.82967cor.test(StrucLome$Leaf, StrucLome$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Leaf and StrucLome$Biomass
## S = 68, p-value = 0.0012
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.81319cor.test(StrucLome$Leaf, StrucLome$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Leaf and StrucLome$Richness
## S = 59.6, p-value = 0.00037
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.83631cor.test(StrucLome$Leaf, StrucLome$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Leaf and StrucLome$Simpson_1.D
## S = 78, p-value = 0.0023
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.78571cor.test(StrucLome$Leaf, StrucLome$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Leaf and StrucLome$BasalHA
## S = 62, p-value = 0.00078
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.82967cor.test(StrucMon$Leaf, StrucMon$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Leaf and StrucMon$TreesHA
## S = 993, p-value = 0.92
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.024781cor.test(StrucMon$Leaf, StrucMon$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Leaf and StrucMon$Biomass
## S = 973, p-value = 0.99
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.0041301cor.test(StrucMon$Leaf, StrucMon$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Leaf and StrucMon$Richness
## S = 908, p-value = 0.8
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.063115cor.test(StrucMon$Leaf, StrucMon$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Leaf and StrucMon$Simpson_1.D
## S = 951, p-value = 0.94
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.01878cor.test(StrucMon$Leaf, StrucMon$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Leaf and StrucMon$BasalHA
## S = 989, p-value = 0.94
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.020704Branch Litter
cor.test(Struc$Branch, Struc$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Branch and Struc$TreesHA
## S = 2907, p-value = 0.021
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.41399cor.test(Struc$Branch, Struc$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Branch and Struc$BasalHA
## S = 2981, p-value = 0.026
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.39903cor.test(Struc$Branch, Struc$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Branch and Struc$Biomass
## S = 3025, p-value = 0.03
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.39016cor.test(Struc$Branch, Struc$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Branch and Struc$Richness
## S = 3359, p-value = 0.077
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.32274cor.test(Struc$Branch, Struc$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Branch and Struc$Simpson_1.D
## S = 3927, p-value = 0.26
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2083Next we will use ANCOVAs to determine if landscapes differ in their relationship between Branch litter and the structural variables.
anova(lm(sqrt(Branch) ~ TreesHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Branch)
## Df Sum Sq Mean Sq F value Pr(>F)
## TreesHA 1 0.165 0.1649 8.20 0.008 **
## Landscape 1 0.006 0.0065 0.32 0.574
## TreesHA:Landscape 1 0.067 0.0670 3.33 0.079 .
## Residuals 27 0.543 0.0201
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Branch) ~ BasalHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Branch)
## Df Sum Sq Mean Sq F value Pr(>F)
## BasalHA 1 0.117 0.1170 5.47 0.027 *
## Landscape 1 0.027 0.0268 1.26 0.272
## BasalHA:Landscape 1 0.060 0.0601 2.81 0.105
## Residuals 27 0.577 0.0214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Branch) ~ Biomass * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Branch)
## Df Sum Sq Mean Sq F value Pr(>F)
## Biomass 1 0.109 0.1085 4.97 0.034 *
## Landscape 1 0.031 0.0310 1.42 0.244
## Biomass:Landscape 1 0.052 0.0518 2.37 0.135
## Residuals 27 0.590 0.0218
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Branch) ~ Richness * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Branch)
## Df Sum Sq Mean Sq F value Pr(>F)
## Richness 1 0.093 0.0925 4.22 0.05 *
## Landscape 1 0.036 0.0356 1.63 0.21
## Richness:Landscape 1 0.061 0.0608 2.77 0.11
## Residuals 27 0.592 0.0219
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Branch) ~ Simpson_1.D * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Branch)
## Df Sum Sq Mean Sq F value Pr(>F)
## Simpson_1.D 1 0.048 0.0484 1.97 0.17
## Landscape 1 0.053 0.0531 2.16 0.15
## Simpson_1.D:Landscape 1 0.015 0.0146 0.59 0.45
## Residuals 27 0.665 0.0246cor.test(StrucLome$Branch, StrucLome$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Branch and StrucLome$TreesHA
## S = 72, p-value = 0.0016
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8022cor.test(StrucLome$Branch, StrucLome$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Branch and StrucLome$Biomass
## S = 56, p-value = 0.00044
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.84615cor.test(StrucLome$Branch, StrucLome$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Branch and StrucLome$Richness
## S = 77.6, p-value = 0.0014
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7868cor.test(StrucLome$Branch, StrucLome$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Branch and StrucLome$Simpson_1.D
## S = 140, p-value = 0.029
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.61538cor.test(StrucLome$Branch, StrucLome$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Branch and StrucLome$BasalHA
## S = 60, p-value = 0.00065
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.83516cor.test(StrucMon$Branch, StrucMon$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Branch and StrucMon$TreesHA
## S = 861, p-value = 0.66
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.11151cor.test(StrucMon$Branch, StrucMon$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Branch and StrucMon$Biomass
## S = 793, p-value = 0.47
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.18172cor.test(StrucMon$Branch, StrucMon$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Branch and StrucMon$Richness
## S = 898, p-value = 0.77
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.073461cor.test(StrucMon$Branch, StrucMon$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Branch and StrucMon$Simpson_1.D
## S = 916, p-value = 0.83
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.054255cor.test(StrucMon$Branch, StrucMon$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Branch and StrucMon$BasalHA
## S = 830, p-value = 0.57
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.14389Detritus Litter
cor.test(Struc$Detritus, Struc$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Detritus and Struc$TreesHA
## S = 3604, p-value = 0.14
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.27344cor.test(Struc$Detritus, Struc$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Detritus and Struc$BasalHA
## S = 3014, p-value = 0.029
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.39237cor.test(Struc$Detritus, Struc$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Detritus and Struc$Biomass
## S = 3035, p-value = 0.031
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.38814cor.test(Struc$Detritus, Struc$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Detritus and Struc$Richness
## S = 3369, p-value = 0.079
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.32072cor.test(Struc$Detritus, Struc$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Detritus and Struc$Simpson_1.D
## S = 3417, p-value = 0.088
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.31114Next we will use ANCOVAs to determine if landscapes differ in their relationship between Detritus litter and the structural variables.
anova(lm(sqrt(Detritus) ~ TreesHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Detritus)
## Df Sum Sq Mean Sq F value Pr(>F)
## TreesHA 1 0.169 0.1690 3.57 0.070 .
## Landscape 1 0.020 0.0197 0.41 0.525
## TreesHA:Landscape 1 0.172 0.1716 3.62 0.068 .
## Residuals 27 1.280 0.0474
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Detritus) ~ BasalHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Detritus)
## Df Sum Sq Mean Sq F value Pr(>F)
## BasalHA 1 0.156 0.1556 3.14 0.088 .
## Landscape 1 0.003 0.0027 0.05 0.816
## BasalHA:Landscape 1 0.145 0.1446 2.92 0.099 .
## Residuals 27 1.337 0.0495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Detritus) ~ Biomass * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Detritus)
## Df Sum Sq Mean Sq F value Pr(>F)
## Biomass 1 0.152 0.1519 3.02 0.094 .
## Landscape 1 0.001 0.0015 0.03 0.866
## Biomass:Landscape 1 0.128 0.1275 2.53 0.123
## Residuals 27 1.359 0.0503
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Detritus) ~ Richness * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Detritus)
## Df Sum Sq Mean Sq F value Pr(>F)
## Richness 1 0.132 0.1321 2.57 0.12
## Landscape 1 0.001 0.0005 0.01 0.92
## Richness:Landscape 1 0.121 0.1213 2.36 0.14
## Residuals 27 1.386 0.0513anova(lm(sqrt(Detritus) ~ Simpson_1.D * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Detritus)
## Df Sum Sq Mean Sq F value Pr(>F)
## Simpson_1.D 1 0.132 0.1318 2.39 0.13
## Landscape 1 0.001 0.0014 0.03 0.88
## Simpson_1.D:Landscape 1 0.017 0.0166 0.30 0.59
## Residuals 27 1.490 0.0552cor.test(StrucLome$Detritus, StrucLome$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Detritus and StrucLome$TreesHA
## S = 82, p-value = 0.0029
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.77473cor.test(StrucLome$Detritus, StrucLome$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Detritus and StrucLome$Biomass
## S = 82, p-value = 0.0029
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.77473cor.test(StrucLome$Detritus, StrucLome$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Detritus and StrucLome$Richness
## S = 78.6, p-value = 0.0015
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.78404cor.test(StrucLome$Detritus, StrucLome$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Detritus and StrucLome$Simpson_1.D
## S = 96, p-value = 0.0058
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.73626cor.test(StrucLome$Detritus, StrucLome$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Detritus and StrucLome$BasalHA
## S = 76, p-value = 0.0021
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.79121cor.test(StrucMon$Detritus, StrucMon$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Detritus and StrucMon$TreesHA
## S = 1053, p-value = 0.73
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.086732cor.test(StrucMon$Detritus, StrucMon$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Detritus and StrucMon$Biomass
## S = 847, p-value = 0.62
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.12597cor.test(StrucMon$Detritus, StrucMon$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Detritus and StrucMon$Richness
## S = 971, p-value = 0.99
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.0020693cor.test(StrucMon$Detritus, StrucMon$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Detritus and StrucMon$Simpson_1.D
## S = 1007, p-value = 0.88
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.039648cor.test(StrucMon$Detritus, StrucMon$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Detritus and StrucMon$BasalHA
## S = 885, p-value = 0.73
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.086957Branch Litter
cor.test(Struc$Branch, Struc$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Branch and Struc$TreesHA
## S = 2907, p-value = 0.021
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.41399cor.test(Struc$Branch, Struc$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Branch and Struc$BasalHA
## S = 2981, p-value = 0.026
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.39903cor.test(Struc$Branch, Struc$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Branch and Struc$Biomass
## S = 3025, p-value = 0.03
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.39016cor.test(Struc$Branch, Struc$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Branch and Struc$Richness
## S = 3359, p-value = 0.077
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.32274cor.test(Struc$Branch, Struc$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Branch and Struc$Simpson_1.D
## S = 3927, p-value = 0.26
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2083Next we will use ANCOVAs to determine if landscapes differ in their relationship between Branch litter and the structural variables.
anova(lm(sqrt(Branch) ~ TreesHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Branch)
## Df Sum Sq Mean Sq F value Pr(>F)
## TreesHA 1 0.165 0.1649 8.20 0.008 **
## Landscape 1 0.006 0.0065 0.32 0.574
## TreesHA:Landscape 1 0.067 0.0670 3.33 0.079 .
## Residuals 27 0.543 0.0201
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Branch) ~ BasalHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Branch)
## Df Sum Sq Mean Sq F value Pr(>F)
## BasalHA 1 0.117 0.1170 5.47 0.027 *
## Landscape 1 0.027 0.0268 1.26 0.272
## BasalHA:Landscape 1 0.060 0.0601 2.81 0.105
## Residuals 27 0.577 0.0214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Branch) ~ Biomass * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Branch)
## Df Sum Sq Mean Sq F value Pr(>F)
## Biomass 1 0.109 0.1085 4.97 0.034 *
## Landscape 1 0.031 0.0310 1.42 0.244
## Biomass:Landscape 1 0.052 0.0518 2.37 0.135
## Residuals 27 0.590 0.0218
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Branch) ~ Richness * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Branch)
## Df Sum Sq Mean Sq F value Pr(>F)
## Richness 1 0.093 0.0925 4.22 0.05 *
## Landscape 1 0.036 0.0356 1.63 0.21
## Richness:Landscape 1 0.061 0.0608 2.77 0.11
## Residuals 27 0.592 0.0219
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Branch) ~ Simpson_1.D * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Branch)
## Df Sum Sq Mean Sq F value Pr(>F)
## Simpson_1.D 1 0.048 0.0484 1.97 0.17
## Landscape 1 0.053 0.0531 2.16 0.15
## Simpson_1.D:Landscape 1 0.015 0.0146 0.59 0.45
## Residuals 27 0.665 0.0246cor.test(StrucLome$Branch, StrucLome$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Branch and StrucLome$TreesHA
## S = 72, p-value = 0.0016
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8022cor.test(StrucLome$Branch, StrucLome$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Branch and StrucLome$Biomass
## S = 56, p-value = 0.00044
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.84615cor.test(StrucLome$Branch, StrucLome$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Branch and StrucLome$Richness
## S = 77.6, p-value = 0.0014
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7868cor.test(StrucLome$Branch, StrucLome$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Branch and StrucLome$Simpson_1.D
## S = 140, p-value = 0.029
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.61538cor.test(StrucLome$Branch, StrucLome$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Branch and StrucLome$BasalHA
## S = 60, p-value = 0.00065
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.83516cor.test(StrucMon$Branch, StrucMon$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Branch and StrucMon$TreesHA
## S = 861, p-value = 0.66
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.11151cor.test(StrucMon$Branch, StrucMon$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Branch and StrucMon$Biomass
## S = 793, p-value = 0.47
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.18172cor.test(StrucMon$Branch, StrucMon$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Branch and StrucMon$Richness
## S = 898, p-value = 0.77
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.073461cor.test(StrucMon$Branch, StrucMon$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Branch and StrucMon$Simpson_1.D
## S = 916, p-value = 0.83
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.054255cor.test(StrucMon$Branch, StrucMon$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Branch and StrucMon$BasalHA
## S = 830, p-value = 0.57
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.14389Flower Litter
cor.test(Struc$Flower, Struc$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Flower and Struc$TreesHA
## S = 3216, p-value = 0.052
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.35153cor.test(Struc$Flower, Struc$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Flower and Struc$BasalHA
## S = 3102, p-value = 0.038
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.37468cor.test(Struc$Flower, Struc$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Flower and Struc$Biomass
## S = 3219, p-value = 0.053
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.35109cor.test(Struc$Flower, Struc$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Flower and Struc$Richness
## S = 2553, p-value = 0.0057
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.48526cor.test(Struc$Flower, Struc$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Flower and Struc$Simpson_1.D
## S = 2573, p-value = 0.0061
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4812Next we will use ANCOVAs to determine if landscapes differ in their relationship between Flower litter and the structural variables.
anova(lm(sqrt(Flower) ~ TreesHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Flower)
## Df Sum Sq Mean Sq F value Pr(>F)
## TreesHA 1 0.0236 0.02359 3.92 0.058 .
## Landscape 1 0.0090 0.00904 1.50 0.231
## TreesHA:Landscape 1 0.0284 0.02840 4.72 0.039 *
## Residuals 27 0.1625 0.00602
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Flower) ~ BasalHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Flower)
## Df Sum Sq Mean Sq F value Pr(>F)
## BasalHA 1 0.0220 0.02201 3.70 0.065 .
## Landscape 1 0.0147 0.01474 2.48 0.127
## BasalHA:Landscape 1 0.0260 0.02602 4.37 0.046 *
## Residuals 27 0.1608 0.00595
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Flower) ~ Biomass * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Flower)
## Df Sum Sq Mean Sq F value Pr(>F)
## Biomass 1 0.0207 0.02068 3.44 0.075 .
## Landscape 1 0.0160 0.01598 2.66 0.115
## Biomass:Landscape 1 0.0246 0.02455 4.08 0.053 .
## Residuals 27 0.1623 0.00601
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Flower) ~ Richness * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Flower)
## Df Sum Sq Mean Sq F value Pr(>F)
## Richness 1 0.0434 0.0434 7.96 0.0089 **
## Landscape 1 0.0164 0.0164 3.00 0.0944 .
## Richness:Landscape 1 0.0166 0.0166 3.04 0.0928 .
## Residuals 27 0.1472 0.0055
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Flower) ~ Simpson_1.D * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Flower)
## Df Sum Sq Mean Sq F value Pr(>F)
## Simpson_1.D 1 0.0238 0.02383 3.70 0.065 .
## Landscape 1 0.0246 0.02461 3.82 0.061 .
## Simpson_1.D:Landscape 1 0.0012 0.00121 0.19 0.669
## Residuals 27 0.1739 0.00644
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.test(StrucLome$Flower, StrucLome$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Flower and StrucLome$TreesHA
## S = 81.2, p-value = 0.0018
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.77686cor.test(StrucLome$Flower, StrucLome$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Flower and StrucLome$Biomass
## S = 73.2, p-value = 0.0011
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7989cor.test(StrucLome$Flower, StrucLome$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Flower and StrucLome$Richness
## S = 58.7, p-value = 0.00034
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.83862cor.test(StrucLome$Flower, StrucLome$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Flower and StrucLome$Simpson_1.D
## S = 91.2, p-value = 0.0032
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.74931cor.test(StrucLome$Flower, StrucLome$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Flower and StrucLome$BasalHA
## S = 81.2, p-value = 0.0018
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.77686cor.test(StrucMon$Flower, StrucMon$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Flower and StrucMon$TreesHA
## S = 977, p-value = 0.97
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.0082687cor.test(StrucMon$Flower, StrucMon$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Flower and StrucMon$Biomass
## S = 969, p-value = 1
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0cor.test(StrucMon$Flower, StrucMon$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Flower and StrucMon$Richness
## S = 695, p-value = 0.26
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.28276cor.test(StrucMon$Flower, StrucMon$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Flower and StrucMon$Simpson_1.D
## S = 587, p-value = 0.11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.39375cor.test(StrucMon$Flower, StrucMon$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Flower and StrucMon$BasalHA
## S = 952, p-value = 0.94
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.017617Fruit Litter
cor.test(Struc$Fruit, Struc$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Fruit and Struc$TreesHA
## S = 5340, p-value = 0.68
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.076636cor.test(Struc$Fruit, Struc$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Fruit and Struc$BasalHA
## S = 4587, p-value = 0.69
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.075255cor.test(Struc$Fruit, Struc$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Fruit and Struc$Biomass
## S = 4495, p-value = 0.62
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.093769cor.test(Struc$Fruit, Struc$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Fruit and Struc$Richness
## S = 4608, p-value = 0.7
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.071054cor.test(Struc$Fruit, Struc$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Fruit and Struc$Simpson_1.D
## S = 4586, p-value = 0.69
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.075369Next we will use ANCOVAs to determine if landscapes differ in their relationship between Fruit litter and the structural variables.
anova(lm(sqrt(Fruit) ~ TreesHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Fruit)
## Df Sum Sq Mean Sq F value Pr(>F)
## TreesHA 1 0.000 0.00027 0.01 0.91
## Landscape 1 0.029 0.02894 1.30 0.26
## TreesHA:Landscape 1 0.001 0.00089 0.04 0.84
## Residuals 27 0.602 0.02230anova(lm(sqrt(Fruit) ~ BasalHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Fruit)
## Df Sum Sq Mean Sq F value Pr(>F)
## BasalHA 1 0.006 0.00554 0.25 0.62
## Landscape 1 0.028 0.02774 1.25 0.27
## BasalHA:Landscape 1 0.001 0.00143 0.06 0.80
## Residuals 27 0.598 0.02213anova(lm(sqrt(Fruit) ~ Biomass * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Fruit)
## Df Sum Sq Mean Sq F value Pr(>F)
## Biomass 1 0.007 0.00707 0.32 0.58
## Landscape 1 0.027 0.02699 1.22 0.28
## Biomass:Landscape 1 0.000 0.00040 0.02 0.89
## Residuals 27 0.598 0.02214anova(lm(sqrt(Fruit) ~ Richness * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Fruit)
## Df Sum Sq Mean Sq F value Pr(>F)
## Richness 1 0.003 0.00339 0.15 0.70
## Landscape 1 0.025 0.02543 1.14 0.29
## Richness:Landscape 1 0.001 0.00139 0.06 0.80
## Residuals 27 0.602 0.02230anova(lm(sqrt(Fruit) ~ Simpson_1.D * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Fruit)
## Df Sum Sq Mean Sq F value Pr(>F)
## Simpson_1.D 1 0.080 0.0803 4.12 0.052 .
## Landscape 1 0.016 0.0157 0.80 0.378
## Simpson_1.D:Landscape 1 0.010 0.0098 0.50 0.485
## Residuals 27 0.526 0.0195
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.test(StrucLome$Fruit, StrucLome$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Fruit and StrucLome$TreesHA
## S = 366, p-value = 0.99
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.0054945cor.test(StrucLome$Fruit, StrucLome$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Fruit and StrucLome$Biomass
## S = 332, p-value = 0.78
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.087912cor.test(StrucLome$Fruit, StrucLome$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Fruit and StrucLome$Richness
## S = 306, p-value = 0.6
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.15956cor.test(StrucLome$Fruit, StrucLome$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Fruit and StrucLome$Simpson_1.D
## S = 348, p-value = 0.89
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.043956cor.test(StrucLome$Fruit, StrucLome$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Fruit and StrucLome$BasalHA
## S = 334, p-value = 0.79
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.082418cor.test(StrucMon$Fruit, StrucMon$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Fruit and StrucMon$TreesHA
## S = 1098, p-value = 0.6
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.1332cor.test(StrucMon$Fruit, StrucMon$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Fruit and StrucMon$Biomass
## S = 911, p-value = 0.81
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.059886cor.test(StrucMon$Fruit, StrucMon$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Fruit and StrucMon$Richness
## S = 982, p-value = 0.96
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.013451cor.test(StrucMon$Fruit, StrucMon$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Fruit and StrucMon$Simpson_1.D
## S = 976, p-value = 0.98
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.0073035cor.test(StrucMon$Fruit, StrucMon$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Fruit and StrucMon$BasalHA
## S = 946, p-value = 0.93
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.02381Total Litter
cor.test(Struc$Total, Struc$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Total and Struc$TreesHA
## S = 4572, p-value = 0.68
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.078242cor.test(Struc$Total, Struc$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Total and Struc$BasalHA
## S = 3535, p-value = 0.12
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.28727cor.test(Struc$Total, Struc$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Total and Struc$Biomass
## S = 3501, p-value = 0.11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.29418cor.test(Struc$Total, Struc$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Total and Struc$Richness
## S = 3585, p-value = 0.13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.27712cor.test(Struc$Total, Struc$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: Struc$Total and Struc$Simpson_1.D
## S = 3392, p-value = 0.083
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.31619Next we will use ANCOVAs to determine if landscapes differ in their relationship between Total litter and the structural variables.
anova(lm(sqrt(Total) ~ TreesHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Total)
## Df Sum Sq Mean Sq F value Pr(>F)
## TreesHA 1 0.275 0.275 2.99 0.09513 .
## Landscape 1 1.375 1.375 14.97 0.00062 ***
## TreesHA:Landscape 1 0.816 0.816 8.89 0.00602 **
## Residuals 27 2.480 0.092
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Total) ~ BasalHA * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Total)
## Df Sum Sq Mean Sq F value Pr(>F)
## BasalHA 1 0.442 0.442 4.35 0.0466 *
## Landscape 1 1.075 1.075 10.58 0.0031 **
## BasalHA:Landscape 1 0.685 0.685 6.74 0.0150 *
## Residuals 27 2.744 0.102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Total) ~ Biomass * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Total)
## Df Sum Sq Mean Sq F value Pr(>F)
## Biomass 1 0.486 0.486 4.63 0.0405 *
## Landscape 1 1.027 1.027 9.79 0.0042 **
## Biomass:Landscape 1 0.602 0.602 5.74 0.0237 *
## Residuals 27 2.831 0.105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Total) ~ Richness * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Total)
## Df Sum Sq Mean Sq F value Pr(>F)
## Richness 1 0.507 0.507 4.67 0.0396 *
## Landscape 1 0.973 0.973 8.97 0.0058 **
## Richness:Landscape 1 0.537 0.537 4.95 0.0346 *
## Residuals 27 2.929 0.108
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(sqrt(Total) ~ Simpson_1.D * Landscape, data = Struc)) ## Analysis of Variance Table
##
## Response: sqrt(Total)
## Df Sum Sq Mean Sq F value Pr(>F)
## Simpson_1.D 1 0.44 0.441 3.23 0.083 .
## Landscape 1 0.77 0.773 5.66 0.025 *
## Simpson_1.D:Landscape 1 0.04 0.045 0.33 0.573
## Residuals 27 3.69 0.137
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.test(StrucLome$Total, StrucLome$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Total and StrucLome$TreesHA
## S = 56, p-value = 0.00044
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.84615cor.test(StrucLome$Total, StrucLome$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Total and StrucLome$Biomass
## S = 58, p-value = 0.00054
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.84066cor.test(StrucLome$Total, StrucLome$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Total and StrucLome$Richness
## S = 52.6, p-value = 0.00019
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.85557cor.test(StrucLome$Total, StrucLome$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Total and StrucLome$Simpson_1.D
## S = 86, p-value = 0.0036
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.76374cor.test(StrucLome$Total, StrucLome$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucLome$Total and StrucLome$BasalHA
## S = 50, p-value = 0.00019
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.86264cor.test(StrucMon$Total, StrucMon$TreesHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Total and StrucMon$TreesHA
## S = 1120, p-value = 0.54
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.15591cor.test(StrucMon$Total, StrucMon$Biomass, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Total and StrucMon$Biomass
## S = 972, p-value = 0.99
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.0030976cor.test(StrucMon$Total, StrucMon$Richness, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Total and StrucMon$Richness
## S = 1041, p-value = 0.77
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.074496cor.test(StrucMon$Total, StrucMon$Simpson_1.D, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Total and StrucMon$Simpson_1.D
## S = 1092, p-value = 0.61
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.12729cor.test(StrucMon$Total, StrucMon$BasalHA, method="spearman")##
## Spearman's rank correlation rho
##
## data: StrucMon$Total and StrucMon$BasalHA
## S = 1020, p-value = 0.84
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.052795Figure 5
Let’s use Ordination techniques to evaluate how plots differentiate.
First we make a copy of Struc and call it
StrucNew StrucNew is then edited to rename the
factor levels of the Succession variable. We then subset
the columns we don’t need, rename StrucNew as
PCAStrucNew, use prcomp to draft the ordination, rename our
landscapes using the list function, relevel using the
relevel function, and graph using
autoplot.
StrucNew <- Struc
levels(StrucNew$Succession) <- list(Early = "DM1", Intermediate = "DM2", `Old-Growth` = "RM",
Early = "DL1", Intermediate = "DL2", `Old-Growth` = "RL")
StrucNew2 <- StrucNew[5:19]
StrucNew3 <- StrucNew2[ , -which(names(StrucNew2) %in% c("Shannon_H","Chao.1", "Family", "Biomass", "Total.Mg.ha.month", "Leaf", "Branch", "Flower", "Fruit", "Detritus"))]
PCAStrucNew <- StrucNew
pca_resx <- prcomp(StrucNew3, scale. = TRUE)
levels(PCAStrucNew$Landscape) <- list(Mountain = "Mountain", Hill = "Lomerio")
PCAStrucNew$Landscape <- relevel(PCAStrucNew$Landscape, ref = "Hill")
autoplot(pca_resx, data = PCAStrucNew, colour = "Landscape", shape = "Succession",
loadings = TRUE, loadings.colour = c('#AA54B4', '#AA54B4', '#AA54B4',
'#AA54B4',
'#B0B315', '#3E99D9'),
loadings.label = TRUE, loadings.label.size = 4, frame = TRUE, frame.type = "norm",
loadings.label.label = c('Tree Density', 'Basal Area', 'Species Density',
'Diversity',
'Total Litter', 'Precipitation'),
loadings.label.colour = c('#AA54B4', '#AA54B4', '#AA54B4',
'#AA54B4',
'#B0B315', '#3E99D9'),
loadings.label.repel=T, loadings.label.hjust = -0.35,
loadings.label.vjust = 0.25) +
guides(fill = "none") +
scale_colour_manual(values=c("darkgreen", "brown")) +
scale_fill_manual(values=c("#23431F", "#67360D")) +
theme_bw(base_size = 15) + guides(color=guide_legend(title="")) +
guides(shape=guide_legend(title="")) +
theme(legend.position = "bottom",
legend.box = "vertical",
legend.spacing.y = unit(0, "pt"),
legend.margin = margin(t = 1, b = 1))
Next, we use the cor function within the eigen function to obtain our Eigvenvalues. We then do some more magic to get in the format we need!
Eigenvalues
eigen(cor(StrucNew3))## eigen() decomposition
## $values
## [1] 3.457685 1.527042 0.520299 0.337626 0.110986 0.046362
##
## $vectors
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -0.497313 0.1774644 0.229850 0.080832 0.750232 0.3145994
## [2,] -0.518787 0.0588992 0.137140 0.267637 -0.096461 -0.7922433
## [3,] -0.514155 0.0058209 0.044145 0.304629 -0.607277 0.5216107
## [4,] -0.408332 -0.1444513 -0.805223 -0.399123 0.063821 -0.0253393
## [5,] -0.208078 -0.6386053 0.488023 -0.552321 -0.075126 -0.0041811
## [6,] 0.095666 -0.7323390 -0.199654 0.603888 0.222138 0.0253084EigenTable2A <- as.data.frame(eigen(cor(StrucNew3))$values)
EigenTable2A <- as.data.frame(EigenTable2A)
colnames(EigenTable2A) <- "Eigenvalue"
EigenTable2A$CumulativePercentVariation <- EigenTable2A$Eigenvalue/sum(EigenTable2A$Eigenvalue)*100
EigenTable2A <- t(EigenTable2A) # "t" function transposes.
EigenTable2B <- as.data.frame(eigen(cor(StrucNew3))$vectors)
EigenTable <- rbind(EigenTable2A,EigenTable2B)
row.names(EigenTable) <-c("Eigenvalue", "Cumulative Percent Variation", colnames(StrucNew3))
colnames(EigenTable) <- c("PC1", "PC2",
"PC3", "PC4", "PC5", "PC6")
kbl(EigenTable) %>% kable_classic()| PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | |
|---|---|---|---|---|---|---|
| Eigenvalue | 3.45769 | 1.52704 | 0.52030 | 0.33763 | 0.11099 | 0.04636 |
| Cumulative Percent Variation | 57.62809 | 25.45070 | 8.67165 | 5.62709 | 1.84977 | 0.77270 |
| TreesHA | -0.49731 | 0.17746 | 0.22985 | 0.08083 | 0.75023 | 0.31460 |
| BasalHA | -0.51879 | 0.05890 | 0.13714 | 0.26764 | -0.09646 | -0.79224 |
| Richness | -0.51416 | 0.00582 | 0.04415 | 0.30463 | -0.60728 | 0.52161 |
| Simpson_1.D | -0.40833 | -0.14445 | -0.80522 | -0.39912 | 0.06382 | -0.02534 |
| Total | -0.20808 | -0.63861 | 0.48802 | -0.55232 | -0.07513 | -0.00418 |
| RainNow | 0.09567 | -0.73234 | -0.19965 | 0.60389 | 0.22214 | 0.02531 |
Structure x Habitat
Do habitats within landscapes differ in structure? We can find out
via anova function on an lm function of each
variable by Age Group.
anova(lm(TreesHA ~ Succession, data = StrucLome))## Analysis of Variance Table
##
## Response: TreesHA
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 344476 172238 48.3 7.3e-06 ***
## Residuals 10 35669 3567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(BasalHA ~ Succession, data = StrucLome))## Analysis of Variance Table
##
## Response: BasalHA
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 865 432 58.6 3e-06 ***
## Residuals 10 74 7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(Simpson_1.D ~ Succession, data = StrucLome))## Analysis of Variance Table
##
## Response: Simpson_1.D
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.139 0.0694 3.19 0.085 .
## Residuals 10 0.218 0.0218
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(Richness ~ Succession, data = StrucLome))## Analysis of Variance Table
##
## Response: Richness
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 8747 4374 62.4 2.2e-06 ***
## Residuals 10 701 70
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(TreesHA ~ Succession, data = StrucMon))## Analysis of Variance Table
##
## Response: TreesHA
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 11916 5958 0.08 0.92
## Residuals 15 1071344 71423anova(lm(BasalHA ~ Succession, data = StrucMon))## Analysis of Variance Table
##
## Response: BasalHA
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 125 62.4 0.52 0.6
## Residuals 15 1789 119.3anova(lm(Simpson_1.D ~ Succession, data = StrucMon))## Analysis of Variance Table
##
## Response: Simpson_1.D
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 0.060 0.0302 0.93 0.42
## Residuals 15 0.488 0.0325anova(lm(Richness ~ Succession, data = StrucMon))## Analysis of Variance Table
##
## Response: Richness
## Df Sum Sq Mean Sq F value Pr(>F)
## Succession 2 1272 636 0.58 0.57
## Residuals 15 16591 1106CONCLUSION: Differences in Hill, not in Mountain! Let’s get some descriptive statistics for the structural variables now.
kbl(summarySE(StrucLome, measurevar="TreesHA", groupvars=c("Succession"),
na.rm = TRUE)[,-4][,-5][,-7], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Succession | N | TreesHA | se |
|---|---|---|---|
| DL1 | 3 | 148.00 | 45.078 |
| DL2 | 6 | 446.67 | 26.385 |
| RL | 4 | 592.00 | 14.697 |
kbl(summarySE(StrucLome, measurevar="BasalHA", groupvars=c("Succession"),
na.rm = TRUE)[,-4][,-5][,-7], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Succession | N | BasalHA | se |
|---|---|---|---|
| DL1 | 3 | 3.3733 | 1.66455 |
| DL2 | 6 | 13.2167 | 1.35706 |
| RL | 4 | 25.5000 | 0.39791 |
kbl(summarySE(StrucLome, measurevar="Simpson_1.D", groupvars=c("Succession"),
na.rm = TRUE)[,-4][,-5][,-7], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Succession | N | Simpson_1.D | se |
|---|---|---|---|
| DL1 | 3 | 0.70800 | 0.18804 |
| DL2 | 6 | 0.93633 | 0.01339 |
| RL | 4 | 0.97200 | 0.00280 |
kbl(summarySE(StrucLome, measurevar="Richness", groupvars=c("Succession"),
na.rm = TRUE)[,-4][,-5][,-7], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Succession | N | Richness | se |
|---|---|---|---|
| DL1 | 3 | 10.667 | 3.1798 |
| DL2 | 6 | 38.667 | 2.9059 |
| RL | 4 | 80.250 | 5.6771 |
kbl(summarySE(StrucMon, measurevar="TreesHA", groupvars=c("Succession"),
na.rm = TRUE)[,-4][,-5][,-7], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Succession | N | TreesHA | se |
|---|---|---|---|
| DM1 | 6 | 547.33 | 132.111 |
| DM2 | 6 | 553.33 | 77.015 |
| RM | 6 | 604.67 | 111.026 |
kbl(summarySE(StrucMon, measurevar="BasalHA", groupvars=c("Succession"),
na.rm = TRUE)[,-4][,-5][,-7], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Succession | N | BasalHA | se |
|---|---|---|---|
| DM1 | 6 | 15.167 | 4.5491 |
| DM2 | 6 | 15.333 | 3.7830 |
| RM | 6 | 20.833 | 4.9626 |
kbl(summarySE(StrucMon, measurevar="Simpson_1.D", groupvars=c("Succession"),
na.rm = TRUE)[,-4][,-5][,-7], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Succession | N | Simpson_1.D | se |
|---|---|---|---|
| DM1 | 6 | 0.82500 | 0.09283 |
| DM2 | 6 | 0.80500 | 0.08586 |
| RM | 6 | 0.93667 | 0.01687 |
kbl(summarySE(StrucMon, measurevar="Richness", groupvars=c("Succession"),
na.rm = TRUE)[,-4][,-5][,-7], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>%
kable_classic()| Succession | N | Richness | se |
|---|---|---|---|
| DM1 | 6 | 42.333 | 13.696 |
| DM2 | 6 | 42.000 | 13.760 |
| RM | 6 | 60.000 | 13.272 |
Species Composition
We can look at which were the dominant species by examineing the
Composition dataframe. First, we will need to join DP and
EF as that’s the classification in the MS. Next, we subset
so that we can obtain Hill and Mountain dataframes.
levels(Composition$Chronosecuence) <- list("EF" = "DP", "EF" = "EF", "IF" = "IF", "OF" = "OF")
SPHill <- droplevels(subset(Composition, Composition$Landscape == "Hill"))
SPMountain <- droplevels(subset(Composition, Composition$Landscape == "Mountain"))Next we can use the aggregate function to sum up the
abundance by Species and Chronosequence. Next
we order the dataframe by abundance.
SPCountH <- as.data.frame(aggregate(abun ~ Species + Chronosecuence, data = SPHill, FUN = sum, na.rm = TRUE))
SPCountH <-SPCountH[order(SPCountH$abun, decreasing = TRUE),]
SPCountH[1:15,]## Species Chronosecuence abun
## 172 Siparuna guianensis EF 242
## 111 Miconia elata EF 160
## 441 Siparuna guianensis IF 119
## 197 Adenocalymma cladotrichum IF 114
## 459 Tapirira guianensis IF 104
## 98 Isertia laevis EF 76
## 2 Adenocalymma cladotrichum EF 74
## 421 Pseudosenefeldera inclinata IF 71
## 80 Henriettea fascicularis EF 69
## 149 Piptocoma discolor EF 65
## 687 Pseudosenefeldera inclinata OF 65
## 306 Henriettea fascicularis IF 59
## 112 Miconia minutiflora EF 55
## 620 Miconia pterocaulon OF 53
## 294 Guatteria punctata IF 49Great, now we can subset by Age Group within hill and take a look.
SPCountH.EF <- droplevels(subset(SPCountH, SPCountH$Chronosecuence == "EF"))
SPCountH.IF <- droplevels(subset(SPCountH, SPCountH$Chronosecuence == "IF"))
SPCountH.OF <- droplevels(subset(SPCountH, SPCountH$Chronosecuence == "OF"))
SPCountH.EF[1:5,]; SPCountH.IF[1:5,]; SPCountH.OF[1:5,]## Species Chronosecuence abun
## 172 Siparuna guianensis EF 242
## 111 Miconia elata EF 160
## 98 Isertia laevis EF 76
## 2 Adenocalymma cladotrichum EF 74
## 80 Henriettea fascicularis EF 69## Species Chronosecuence abun
## 441 Siparuna guianensis IF 119
## 197 Adenocalymma cladotrichum IF 114
## 459 Tapirira guianensis IF 104
## 421 Pseudosenefeldera inclinata IF 71
## 306 Henriettea fascicularis IF 59## Species Chronosecuence abun
## 687 Pseudosenefeldera inclinata OF 65
## 620 Miconia pterocaulon OF 53
## 729 Virola elongata OF 40
## 623 Miconia splendens OF 37
## 649 Oenocarpus bataua OF 27And now we do the same for Mountain.
SPCountM <- as.data.frame(aggregate(abun ~ Species + Chronosecuence,
data = SPMountain, FUN = sum, na.rm = TRUE))
SPCountM <-SPCountM[order(SPCountM$abun, decreasing = TRUE),]
SPCountM[1:15,]## Species Chronosecuence abun
## 64 Henriettea fascicularis EF 459
## 936 Pseudosenefeldera inclinata OF 168
## 336 Graffenrieda conostegioides IF 163
## 1020 Wettinia praemorsa OF 130
## 106 Miconia minutiflora EF 114
## 104 Miconia lourtegiana EF 103
## 1014 Virola elongata OF 73
## 47 Eugenia lambertiana EF 68
## 80 Inga thibaudiana EF 66
## 187 Vismia baccifera EF 66
## 240 Casearia arborea IF 62
## 121 Myrcia splenden EF 60
## 795 Ladenbergia muzonensis OF 56
## 16 Bellucia grossularioides EF 52
## 736 Graffenrieda colombiana OF 50SPCountM.EF <- droplevels(subset(SPCountM, SPCountM$Chronosecuence == "EF"))
SPCountM.IF <- droplevels(subset(SPCountM, SPCountM$Chronosecuence == "IF"))
SPCountM.OF <- droplevels(subset(SPCountM, SPCountM$Chronosecuence == "OF"))
SPCountM.EF[1:5,]; SPCountM.IF[1:5,]; SPCountM.OF[1:5,]## Species Chronosecuence abun
## 64 Henriettea fascicularis EF 459
## 106 Miconia minutiflora EF 114
## 104 Miconia lourtegiana EF 103
## 47 Eugenia lambertiana EF 68
## 80 Inga thibaudiana EF 66## Species Chronosecuence abun
## 336 Graffenrieda conostegioides IF 163
## 240 Casearia arborea IF 62
## 418 Matayba inelegans IF 46
## 279 Cyathea lasiosora IF 45
## 394 Ladenbergia oblongifolia IF 41## Species Chronosecuence abun
## 936 Pseudosenefeldera inclinata OF 168
## 1020 Wettinia praemorsa OF 130
## 1014 Virola elongata OF 73
## 795 Ladenbergia muzonensis OF 56
## 736 Graffenrieda colombiana OF 50Regressions
Let’s calculate regressions for Stand Age vs. Leaf Litter Production. First we need to manually insert the ages of the plots. These must be inserted in this order.
Struc$Age <- as.numeric(c("12", "7", "10", "25", "35", "30", "28", "30", "30",
"12", "5", "5", "17", "5", "13", "26", "35", "31", "37", "22",
"24", "40", "57", "49", "55", "47", "55", "60", "65", "68", "49"))Next we relevel and rename.
levels(Struc$Landscape) <- list(Hill = "Lomerio", Mountain = "Mountain")
Struc$`Stand Age Group` <- Struc$Succession
levels(Struc$`Stand Age Group`) <- list(EH = "DL1", IH = "DL2", OH = "RL",
EM = "DM1", IM = "DM2", OM = "RM")Figure 3
We’re now ready for graphing! For this we will use
ggplot and ggarrange.
GGNew1 <- ggplot(Struc, aes(x=Age, y=Leaf, group = Landscape, color = Landscape)) +
geom_point(size=1, alpha=0.9) +
geom_smooth(size = 0.5, method=lm, aes(fill = Landscape), color = "black") +
scale_color_manual(values=c("#82C27B", "#FBBA84")) +
scale_fill_manual(values=c("#82C27B", "#FBBA84")) +
ylab(bquote('Leaf Litter Production '(Mg~ha^-1~mth^-1))) +
xlab("Stand Age") +
theme_bw(base_size = 9) +
guides(color = "none") +
geom_vline(xintercept = c(20, 40), linetype = "dashed", color = "black")
GGNew2 <- ggplot(Struc, aes(x=Age, y=Total, group = Landscape , color = Landscape)) +
geom_point(size=1, alpha=0.9) +
geom_smooth(size = 0.5, method=lm, aes(fill = Landscape), color = "black") +
scale_color_manual(values=c("#82C27B", "#FBBA84")) +
scale_fill_manual(values=c("#82C27B", "#FBBA84")) +
ylab(bquote('Total Litter Production '(Mg~ha^-1~mth^-1))) +
xlab("Stand Age") +
theme_bw(base_size = 9) +
guides(color = "none") +
geom_vline(xintercept = c(20, 40), linetype = "dashed", color = "black")
ggarrange(GGNew2, GGNew1,
align='h', labels = "auto", legend = "right",
common.legend = T, nrow = 2, ncol = 1)
Session Report
## Analyses were conducted using the R Statistical language (version 4.2.2; R Core
## Team, 2022) on macOS Ventura 13.2.1, using the packages ggthemes (version
## 4.2.4; Arnold J, 2021), lme4 (version 1.1.31; Bates D et al., 2015), Matrix
## (version 1.5.1; Bates D et al., 2022), conover.test (version 1.1.5; Dinno A,
## 2017), lubridate (version 1.9.2; Grolemund G, Wickham H, 2011), Rmisc (version
## 1.5.1; Hope RM, 2022), ggpubr (version 0.6.0; Kassambara A, 2023), lmerTest
## (version 3.1.3; Kuznetsova A et al., 2017), emmeans (version 1.8.4.1; Lenth R,
## 2023), report (version 0.5.6; Makowski D et al., 2023), patchwork (version
## 1.1.2; Pedersen T, 2022), lattice (version 0.20.45; Sarkar D, 2008), ggfortify
## (version 0.4.15; Tang Y et al., 2016), reshape (version 0.8.9; Wickham H,
## 2007), plyr (version 1.8.8; Wickham H, 2011), ggplot2 (version 3.4.1; Wickham
## H, 2016), dplyr (version 1.1.0; Wickham H et al., 2023), zoo (version 1.8.11;
## Zeileis A, Grothendieck G, 2005), lmtest (version 0.9.40; Zeileis A, Hothorn T,
## 2002) and kableExtra (version 1.3.4; Zhu H, 2021).
##
## References
## ----------
## - Arnold J (2021). _ggthemes: Extra Themes, Scales and Geoms for 'ggplot2'_. R
## package version 4.2.4, <https://CRAN.R-project.org/package=ggthemes>.
## - Bates D, Mächler M, Bolker B, Walker S (2015). "Fitting Linear Mixed-Effects
## Models Using lme4." _Journal of Statistical Software_, *67*(1), 1-48.
## doi:10.18637/jss.v067.i01 <https://doi.org/10.18637/jss.v067.i01>.
## - Bates D, Maechler M, Jagan M (2022). _Matrix: Sparse and Dense Matrix Classes
## and Methods_. R package version 1.5-1,
## <https://CRAN.R-project.org/package=Matrix>.
## - Dinno A (2017). _conover.test: Conover-Iman Test of Multiple Comparisons
## Using Rank Sums_. R package version 1.1.5,
## <https://CRAN.R-project.org/package=conover.test>.
## - Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate."
## _Journal of Statistical Software_, *40*(3), 1-25.
## <https://www.jstatsoft.org/v40/i03/>.
## - Hope RM (2022). _Rmisc: Ryan Miscellaneous_. R package version 1.5.1,
## <https://CRAN.R-project.org/package=Rmisc>.
## - Kassambara A (2023). _ggpubr: 'ggplot2' Based Publication Ready Plots_. R
## package version 0.6.0, <https://CRAN.R-project.org/package=ggpubr>.
## - Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest Package: Tests
## in Linear Mixed Effects Models." _Journal of Statistical Software_, *82*(13),
## 1-26. doi:10.18637/jss.v082.i13 <https://doi.org/10.18637/jss.v082.i13>.
## - Lenth R (2023). _emmeans: Estimated Marginal Means, aka Least-Squares Means_.
## R package version 1.8.4-1, <https://CRAN.R-project.org/package=emmeans>.
## - Makowski D, Lüdecke D, Patil I, Thériault R (2023). "Automated Results
## Reporting as a Practical Tool to Improve Reproducibility and Methodological
## Best Practices Adoption." _CRAN_. <https://easystats.github.io/report/>.
## - Pedersen T (2022). _patchwork: The Composer of Plots_. R package version
## 1.1.2, <https://CRAN.R-project.org/package=patchwork>.
## - R Core Team (2022). _R: A Language and Environment for Statistical
## Computing_. R Foundation for Statistical Computing, Vienna, Austria.
## <https://www.R-project.org/>.
## - Sarkar D (2008). _Lattice: Multivariate Data Visualization with R_. Springer,
## New York. ISBN 978-0-387-75968-5, <http://lmdvr.r-forge.r-project.org>.
## - Tang Y, Horikoshi M, Li W (2016). "ggfortify: Unified Interface to Visualize
## Statistical Result of Popular R Packages." _The R Journal_, *8*(2), 474-485.
## doi:10.32614/RJ-2016-060 <https://doi.org/10.32614/RJ-2016-060>,
## <https://doi.org/10.32614/RJ-2016-060>. Horikoshi M, Tang Y (2018). _ggfortify:
## Data Visualization Tools for Statistical Analysis Results_.
## <https://CRAN.R-project.org/package=ggfortify>.
## - Wickham H (2007). "Reshaping data with the reshape package." _Journal of
## Statistical Software_, *21*(12). <https://www.jstatsoft.org/v21/i12/>.
## - Wickham H (2011). "The Split-Apply-Combine Strategy for Data Analysis."
## _Journal of Statistical Software_, *40*(1), 1-29.
## <https://www.jstatsoft.org/v40/i01/>.
## - Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_.
## Springer-Verlag New York. ISBN 978-3-319-24277-4,
## <https://ggplot2.tidyverse.org>.
## - Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar
## of Data Manipulation_. R package version 1.1.0,
## <https://CRAN.R-project.org/package=dplyr>.
## - Zeileis A, Grothendieck G (2005). "zoo: S3 Infrastructure for Regular and
## Irregular Time Series." _Journal of Statistical Software_, *14*(6), 1-27.
## doi:10.18637/jss.v014.i06 <https://doi.org/10.18637/jss.v014.i06>.
## - Zeileis A, Hothorn T (2002). "Diagnostic Checking in Regression
## Relationships." _R News_, *2*(3), 7-10.
## <https://CRAN.R-project.org/doc/Rnews/>.
## - Zhu H (2021). _kableExtra: Construct Complex Table with 'kable' and Pipe
## Syntax_. R package version 1.3.4,
## <https://CRAN.R-project.org/package=kableExtra>.